Abstract.A very simple and accurate numerical method which is applicable to systems of differentio-integral equations with quite general boundary conditions has been devised. Although the basic idea of this method stems from the Keller Box method, it solves the problem of systems of differential equations involving integral operators not previously considered by the Keller Box method. Two main preparatory stages are required: (i) a merging procedure for differential equations and conditions without integral operators and; (ii) a reduction procedure for differential equations and conditions with integral operators. The differencing processes are effectively simplified by means of the unit-step function. The nonlinear difference equations are solved by Newton's method using an efficient block arrow-like matrix factorization technique. As an example of the application of this method, the systems of equations for combined gravity body force and forced convection in laminar film condensation can be solved for prescribed values of physical constants.