In one numerical method for integration of the diffusion equation in one dimension, the time derivative is replaced by a finite difference in a time interval, and the space derivative by the mean of its values at the beginning and end of the interval. This leads to a set of ordinary differential equations, one for each interval, which have to be solved in succession. Each of these equations is second-order with two-point boundary conditions; the process of integration from one end is severely unstable, the more so the smaller the tune interval. This paper is concerned with a practical, direct and stable method for solving them by integration of two first-order equations, one being integrated inwards and the other outwards, one boundary condition being satisfied in each integration. The extension to axially symmetrical diffusion is briefly considered.