When you add trust-region modifications to Newton-Krylov, you have to constrain the trust-region optimization problem to the Krylov subspace. It’s just a bit of extra low-d linear algebra, but it ends up pushing the need for knowledge of the Krylov subspace into the trust-region calculation. You can’t just stick a trust-region algorithm and a Krylov linear algebra solver together and have them “just work.” So I expect a library like NLSolve will need some modification to work with Krylov methods. Unless of course NLSolve has already has that Krylov tweak incorporated.
And revelt’s PseudoArclengthContinuation.jl has already put these algorithms together.