Consider a nonlinear operator equation in the form: F(x)=0 (1.1) on a pair of Hilbert spaces. Assume that equation (1.1) has a solution, not necessarily unique, and the operator F is twice Fr\'echet differentiable without such structural assumptions as monotonicity, invertibility of F'(x) etc. To avoid the ill-posed inversion of the Fr\'echet derivative operator F'(x) various discrete and continuous methods based on a regularization were suggested. A principal point in the numerical implementation of regularized Newton's and Gauss-Newton's procedures is the computation of the operators (F'(x)+\alpha I)^{-1} and (F^{\prime*}(x)F'(x)+\alpha I)^{-1} respectively. This computation for certain operators requires a considerable effort in many applications. Besides it may decrease the accuracy of the approximate solution. In order to deal with it, a novel iteratively regularized algorithm with simultaneous updates of the operator (F^{\prime*}(x_n)F'(x_n)+\alpha_n I)^{-1} is proposed: x_{n+1}=x_n-B_n [F^{\prime*}(x_n)F(x_n)+\alpha_n(x_n-x_0)], B_{n+1}=[I-\lambda (F^{\prime*}(x_n)F'(x_n)+\alpha_n I)]B_n+\lambda I, A convergence theorem is proved. The stability of the process towards noise in the data is analyzed, and a stopping time is chosen so that the method converges as the noise level tends to zero. The scheme is illustrated by a numerical example in which a nonlinear inverse problem of gravitational sounding is investigated. Based on theoretical and numerical results the recommendations on the choice of \alpha_n, \lambda and B_0 are given. ~