As per the documentation of scipy newton_krylov method, the method can be updated any scipy.sparse.linalg linear solver

Krylov method to use to approximate the Jacobian. Can be a string, or a function implementing the same interface as the iterative solvers in scipy.sparse.linalg. If a string, needs to be one of: ‘lgmres’, ‘gmres’, ‘bicgstab’, ‘cgs’, ‘minres’, ‘tfqmr’.

I tried updating the method to spsolve

from scipy.sparse.linalg import spsolve
r = newton_krylov(Res_H, xin=H, method=spsolve)

but I was getting an error message as below

~Anaconda3libsite-packagesscipyoptimizenonlin.py in nonlin_solve(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)
    307         # The tolerance, as computed for scipy.sparse.linalg.* routines
    308         tol = min(eta, eta*Fx_norm)
--> 309         dx = -jacobian.solve(Fx, tol=tol)
    310 
    311         if norm(dx) == 0:

~Anaconda3libsite-packagesscipyoptimizenonlin.py in solve(self, rhs, tol)
   1574             sol, info = self.method(self.op, rhs, **self.method_kw)
   1575         else:
-> 1576             sol, info = self.method(self.op, rhs, tol=tol, **self.method_kw)
   1577         return sol
   1578 

Just thinking out loud, Can this actually be done WITHOUT modifying the core code ?
If so can someone please shed some light ?