Monotone iteration procedures are developed for the solution of a class of coupled and highly nonlinear boundary value problems. These procedures are of guaranteed convergence and can be used for a numerical investigation of the multiplicity structure of steady-state solutions. The monotone iterations provide also a basis for an automatic scheme for adapting and refining the collocation point grid between iterations. The local error of approximation is controlled very accurately and problems exhibiting extremely steep multiple concentration profiles of significantly different scales can be treated.