Abstract In this work, we propose a free derivative quasi-Newton method for solving large-scale nonlinear systems of equation. We introduce a two-stage linear search direction and develop its global convergence theory. Besides, we prove that the method enjoys superlinear convergence rate. Finally, numerical experiments illustrate that the proposed method is competitive with respect to Newton-Krylov methods and other well-known methods for solving large-scale nonlinear systems of equations.