Abstract: The numerical simulation based on the mixed form of Richards’ equation has the mass conservative property, and can be used to simulate the variably saturated seepage which is prior to other forms of equations. However, on material interface, the distribution of moisture content which is required for the FEM numerical simulation is discontinuous, thus the mixed form of Richards’ equation cannot be utilized in the FEM simulation for heterogeneous soil. In this research, a numerical method was proposed to simulate the variably saturated water flow in heterogeneous soil based on the mixed form of the equation. In the proposed method, two nodal numbers are assigned to each node which locates on the material interface so that the moisture content of each node can be expanded in a truncated Taylor series with respect to pressure head with different material properties. The continuity of the total head is imposed by the Lagrange multiplier method. The proposed method and program based on Matlab were validated against two examples stemmed from other literature. One example was a 1 d problem for the infiltration and draining process of a two layered soil column, and another example was a 2 d infiltration problem for a two layered soil column. And the verification indicated a good agreement between our modeling results and analytical solution for 1 d problem as well as the numerical results for 2 d problem.