Fixed point iterations \(u^{(\nu + 1)}(\lambda) = F(u^{(\nu)},\lambda)\), \(\nu = 0,1,2,\dots,\) may fail to converge against a fixed point \(u^*(\lambda) = F(u^*(\lambda),\lambda)\) of the parameter dependent nonlinear \(F: \mathbb{R}^ N \times \mathbb{R} \to \mathbb{R}^ N\) due to a few unstable eigenvalues of \(F^*_ u := F_ u(u^*(\lambda),\lambda)\) outside the unit circle. A stabilized iteration is proposed which is a combination of Newton's method on the unstable eigenspace (which is assumed to be of small dimension) and a fixed point iteration on its orthogonal complement. The use is recommended in the context of continuation where the number of unstable eigenvalues varies with the parameter \(\lambda\). A similar idea has been used by \textit{H. Jarausch} and \textit{W. Mackens} [Numer. Math. 50, 633-653 (1987; Zbl 0647.65036)] for symmetric \(F_ u\) and is extended here to nonsymmetric Jacobians. A convergence analysis for this stabilized iteration is performed. Details of an algorithm called recursive projection method are given describing the approximation and updating of the invariant unstable eigenspaces during the continuation process. For large dimensional systems arising in discretizations of partial differential equations this algorithm is reported to be highly efficient in comparison with classical pseudo-arclength algorithms. Numerical experiments are performed for the computation of stable and unstable steady state solutions of evolution equations \(u_ t = G(u,\lambda)\) where \(F\) is given by explicit Runge-Kutta steps with \(N = 40\) spatial grid points.