An efficient recursive procedure to solve N-dimensional nonlinear equations using the modified Broyden method is described, This procedure is extended to include the direct inversion in the iterative subspace (DIIS) method to further improve the rate of convergence in iterative calculations. In the recursive procedures. the approximate solutions are constructed as linear combinations of n vectors of length N. The calculations are reduced to determine the appropriate coefficients of the linear combinations. The coefficients are evaluated through small matrix operations, the size of which are at most (n + 1) x (n + 1) except for the generation of a n x n matrix, where n is the iteration number. Storage is required only for the n vectors and the small matrices. The procedures described below can be applied to large systems. To examine the efficiency of the methods, some numerical results are presented in the context of self-consistent calculations of liquid structure using the reference interaction site model (RISM) integral equation and a molecule-site form of the Ornstein-Zernike integral equation. The results indicate that significant acceleration with respect to the Picard iteration method has been achieved by the recursive procedures: The converged solution is obtained in a very small number of iterations and in a fraction of the CPU time. Moreover, the extended method which includes the DIIS approach has further improved the rate of convergence. It also enables solutions to be obtained in an otherwise divergent case.