The multi-aquifer system separated by semipervious leaky beds was analyzed. The finite difference scheme of the Crank-Nicolson method is applied to obtain the solution for this system. The solution of this scheme was compared with the analytical solution for two-layer aquifer systems with one-dimensional steady state. The results showed a good agreement between analytical and numerical solution for two-layer aquifer systems. So, the numerical scheme can be extended to multi-aquifer system. When the pumping is tried for single or multi aquifer, the computation of the groundwater heads was possible for each aquifer in the multi-aquifer with two-dimensional system. So, it might be helpful for the effective groundwater management.