I would guess that the reason why adding an L1 penalty slows things down significantly is that an L1 penalty is not differentiable (i.e. absolute value), while the L2 penalty is. This means that the surface of the function will not be smooth, and so standard quasi-Newton's methods will have a lot of trouble with this problems. Recall that one way to think of a quasi-Newton's method is that it makes a quadratic approximation of the function and then the initial proposal will the maximum of that approximation. If the quadratic approximation matches fairly well to the target function, we should expect the proposal to be close the maximum (or minimum, depending on how you look at the world). But if your target function is non-differentialable, this quadratic approximation may be very bad, thus taking many more iterations to converge (assuming it does).
If you've found an R-package that implements BFGS for L1 penalties, by all means try it. BFGS, in general, is a very generic algorithm for optimization. As is the case with any generic algorithm, there will be plenty of special cases where it does not do well. Algorithms that are specially tailored to your problem clearly should do better (assuming the package is as good as it's author claims: I haven't heard of lbfgs, but there's a whole lot of great things I haven't heard of. Update: I have used R's lbfgs package, and the L-BFGS implementation they have is quite good! I still haven't used it's OWL-QN algorithm, which is what the OP is referring to).
If it doesn't work out for you, you might want to try the "Nelder-Mead" method with R's optim. It does not use derivatives for optimization. As such, it will typically be slower on a smooth function but stabler on an unsmooth function such as you have.