Effectively simulating unsaturated flow is of great significance in many areas such as the soil slope stability analysis, seepage of earth-rockfill dam, contaminant transport. Due to the strongly nonlinear characteristics of Richards equation for the unsaturated flow, numerical solution schemes such as the finite element method with an efficient iterative algorithm usually have to be employed. Picard method as a practical nonlinear iterative method is widely applied to the unsaturated flow field; but usually suffers from convergence oscillation, slow convergence rate and inaccurate solution. In order to improve the computing performance, an efficient finite element procedure based on the adaptive relaxed Picard method is developed. Through the simulations of 1D and 2D unsaturated seepage problems, accuracy, efficiency and robustness of the proposed procedure are validated by comparing with the traditional methods. It is shown that the adaptive relaxed Picard method can effectively reduce the convergence oscillation and significantly improve the convergence rate with the accuracy guaranteed. The proposed procedure provides a helpful reference for the program development and application to unsaturated flow.