According to the STAN homepage, STAN is capable of penalized maximum likelihood (BFGS) optimization. I am using R package rstan but I haven't found any way how to use this method. I tried to look at the ?stan help for the stan() function, but the only available options algorithms are "NUTS" and "HMC".
I am using rstan version 2.5.0.
You want to look at
help(optimizing, package = "rstan"), i.e. theoptimizing()function in the RStan package.