BFGS is not a nonlinear solver method. It seeks to maintain a SPD Hessian for optimization. The thing to remember is that optimization and nonlinear solvers are not the same thing. Using an algorithm desinged for one of them to solve the other is a bad, bad, bad idea.
Broyden looks for a general low-rank change to an approximate Jacobian and is designed for nonlinear solver applications. Broyden needs a good preconditioner to perform well and has fallen out of favor in the last few decades. Newton-Krylov methods do as well and are more common in modern solvers. The central issue is that quasi-Newton methods, when implemented correctly, consume storage during the nonlinear iteration. Newton-GMRES, for example, consumes storage in the linear iteration and that is far easier to manage and makes nonlinear convergence a simpler problem.
My package (SIAMFANLEquations.jl) has Newton-GMRES working, but the documentation and examples are not really ready for prime time. I’ll have something decent in a month and everything in the master branch works correctly. The latest released version (.0.3.0) works, but the memory allocations are better by far in the master branch.