In this work new quasi-Newton methods for solving large-scale nonlinear systems of equations are presented. In these methods q(>1) columns of the approximation of the inverse Jacobian matrix are updated in such a way that the q last secant equations are satisfied (whenever possible) at every iteration. An optimal maximum value for q that makes the method competitive is strongly suggested. The best implementation from the point of view of linear algebra and numerical stability is proposed and a local convergence result for the case q=2 is proved. Several numerical comparative tests with other quasi-Newton methods are carried out.