Abstract: There are often some zero coefficients in the coefficient matrix after the multi-dimensional saturated-unsaturated flow equation is discretized, which will affect the storage space and calculation efficiency in the solution process. A reasonable coefficient matrix optimization method can not only guarantee the simulation accuracy, but also improve the calculation efficiency and universality of the model when it is applied in large-scale grid cells. This paper proposed a generalized multi-dimensional saturated-unsaturated flow model based on spatial linkers. Finite difference method is used to solve the multi-dimensional saturated-unsaturated flow equation. The connections between each two effective units in different directions were pre-recorded, which will avoid the calculation and storage of zero coefficients. The matrix marking method is applied to further simplify the coefficient matrix. The simulation accuracy and efficiency of the model are verified by simulating four classic cases, such as groundwater level rise and seepage surface drainage, as well as three long periods of light rain infiltration in the field. The comparison results show that the simulation accuracy of the model is comparable to some mature softwares, such as Hydrus and VSF, under multiple dimensions and different boundary conditions (including natural rainfall, seepage surface, etc.), while the efficiency is slightly lower than that with Hydrus. Field simulation results also show that the parameter n in VG model can make greatest changes in simulation results with smallest variation and needs to be calibrated in priority. Because the effect of macropore flow has not been considered in the model, the soil moisture response time of each layer lags behind the measured process. At the end of the three rainfall calculation periods, more than 80% of the infiltration water is still stored in the upper 40 cm soil layer, and only about 4%—12% of the infiltration water has been converted to phreatic water. The model in this paper is expected to be an important supplement to the traditional multi-dimensional saturated-unsaturated flow models.