technique. If a square net (h X h) is constructed over R + F, the coefficients of a pair of triangular matrices can be so chosen that the product of the matrices, operating on the values of U at the nodes of the net, approximates to the differential equation (or boundary condition) with truncation error of order 0(h) at every node. Indeed, the coefficients of the triangular matrices may be evaluated at each internal node in turn by solving a set of six linear algebraic equations for each node. Once these coefficients have been computed, the differential equation may be solved numerically for any set of boundary values (and for any function G(x, y)) by back-substitution of the pair of triangular matrices. For any particular set of boundary values, this double back-substitution requires altogether only six multiplications per node.