A numerical solver developed for the solution of parabolic partial differential equations involving two spatial scales is presented. The equations are discretized using the finite volume method, and the resulting system of ordinary differential equations is solved using the stiff code DASSL. The solver has been implemented in the form of a collection of FORTRAN subroutines. The large sparse linear systems which occur during the solution process are solved by a direct efficient method based on a complete LU factorization. The method is fully described in the paper. The storage space needed by the algorithm is reduced by a factor of 3m/4, m being the number of micro grid points, with respect to the standard LINPACK band solver. A comparison of the number of floating point operations involved in the different algorithms shows that this approach reduces the operation count for the factorization phase by a factor of about m, with respect to the LINPACK solver. For the substitution phase the improvement ratio is approximately equal to 3m/4.