Is there a julia package that uses Newton-Krylov method for systems of nonlinear equations?

Yes, I’ve mentioned this to @pkofod a few times.

I guess it is mostly there. You can pass the linearsolver with the argument linsolve