The local rectangular refinement (LRR) solution-adaptive gridding method automatically produces orthogonal unstructured adaptive grids and incorporates multiple-scale finite differences to discretize systems of elliptic governing partial differential equations (PDEs). The coupled non-linear discretized equations are solved simultaneously via Newton's method with a Bi-CGSTAB linear system solver. The grids' unstructured nature produces a nonstandard sparsity pattern within the Jacobian. The effects of two discretization schemes (LRR multiple-scale stencils and traditional single-scale stencils) on Jacobian bandwidth, convergence speed, and solution accuracy are studied. With various point orderings, for two simple problems with analytical solutions, the LRR multiple-scale stencils are seen to: (1) produce Jacobians of smaller bandwidths than those resulting from the traditional single-scale stencils; (2) lead to significantly faster Newton's method convergence than the single-scale stencils; and (3) produce more accurate solutions than the single-scale stencils. The LRR method, including the LRR multiple-scale stencils, is finally applied to an engineering problem governed by strongly coupled, highly non-linear PDEs: a steady-state lean Bunsen flame with complex chemistry, multicomponent transport, and radiation modeling. Very good agreement is observed between the computed flame height and previously published experimental data. Copyright (C) 2001 John Wiley & Sons, Ltd.