The unsaturated zone plays an important role in groundwater recharge, discharge, and the ecological environment. Mathematical models are essential to the investigations of flow processes in the unsaturated zone. The governing equations of the mathematic models for unsaturated flow are nonlinear since their unsaturated hydraulic parameters of the equations depend on their solutions (i.e., pressure or moisture content). The nonlinearity can become very strong with flow in relatively dry soils, and the equation has hyperbolic characteristics. Therefore, the traditional numerical methods, such as the finite difference method and the finite element method, could lead to the numerical oscillation, dispersion, and the divergence of the solution if improper time and space steps were selected. In this study, a finite analytic method (FAM) to solve Richards’ equation was developed. The basic idea of the FAM is the incorporation of the analytic solution in a small local element to formulate the algebraic representation of the partial differential equation of the unsaturated flow. Therefore, the finite analytic numerical method can effectively control the numerical oscillation and dispersion. The convergence and stability of finite analytic numerical scheme are then proven by a rigorous mathematical analysis. Numerical experiments are then used to show that FAM is highly accurate by comparing its results with those obtained using analytical solutions by a modified Picard finite difference method. In addition, we also demonstrate that FAM reproduces results of laboratory experiments. Therefore, it can be considered an appropriate simulation method for the unsaturated flow.