A numerical procedure is described for the computation of a two-layer horizontal groundwater flow in a homogenous aquifer overlying a horizontal impermeable base. The procedure can be extended to more than two layers and to any base configuration. The fluid density is assumed to be constant thoughout each layer. The numerical procedure has second-order accuracy and is unconditionally stable. As it contains no iterative routines, it is also exceptionally economical in computing time and memory requirements.