In the analysis of boundary layers by the local nonsimilarity solution method, the central task is the numerical solution of a set of simultaneous ordinary differential equations. The current method of solving these equations is forward integration in conjunction with a shooting method for determining certain of the starting values for the integration. This approach has been found to be less and less effective as the number of simultaneous equations increases. The numerical solution scheme described here is able to deal effectively with the multiequation systems encountered in local nonsimilarity boundary-layer analysis. It employs integrated forms of the governing differential equations. The key feature of the integrated forms is that they already satisfy the boundary conditions. With these equations, initial, almost arbitrary guesses of the profiles are refined automatically until convergence is attained. For concreteness, the description of the scheme is illustrated by a nonsimilar natural convection boundary-layer problem.