The computation of solution paths for continuation problems requires the solution of a sequence of nonlinear systems of equations. Each nonlinear system can be solved by computing the solution of a succession of linear systems of equations determined by Jacobian matrices associated with the nonlinear system of equations. Points on the solution path where the Jacobian matrix is singular are referred to as singular points and require special handling. They may be turning points or bifurcation points. In order to detect singular points, it is essential to monitor the eigenvalues of smallest magnitude of the Jacobian matrices generated as the solution path is traversed. We describe iterative methods for the computation of solution paths for continuation problems so large that factorization of the Jacobian matrices is infeasible or impractical. The iterative methods simultaneously solve linear systems of equations determined by the Jacobian matrices and compute a few eigenvalue–eigenvector pairs associated with the eigenvalues of smallest magnitude of each Jacobian. A bordering algorithm with a pseudo-arclength parametrization is applied in the vicinity of turning points to overcome the singularity of the Jacobian. A bifurcation perturbation strategy is used to compute solution paths at bifurcation points. Our iterative methods are based on the block-Lanczos algorithm and are applicable to problems with large symmetric Jacobian matrices.