Variably saturated flow in an unsaturated zone is a critical subject due to its diverse applications in earth science, waste management, agriculture, and geotechnical problems. The governing Richards equation (RE) of the mathematical models for unsaturated flow is very nonlinear, attracting researchers to solve it quickly and accurately. In this paper, a numerical model is proposed to solve mixed form of RE. In time, backward Euler method is implemented to discrete the equation, and in space, the equation is discretized by the finite difference method. Because of the dependence of both water content and unsaturated hydraulic conductivity on soil water pressure head, the equation must be solved iteratively using the Picard scheme. It is assumed that initial suction head is 1000 cm and suction head at top boundary is 75 cm throughout the simulation. The model is first tested for a homogenous soil water system using numerical simulation data from the literature. The comparison between simulated and observed value demonstrate that the new algorithm is robust and practical with fast convergence. Later, model is extended for layered soil water system assuming a constant pressure head surface boundary condition and van Genuchten parameter is used to define the soil constitutive relationship for each soil layer. The model shows promising results, allowing it to be used for further applications, such as managed aquifer recharge, plant water uptake, and coupled with the advection-dispersion transport equation to predict solute movement.