Abstract: To solve complicated mechanics problems in geotechnical engineering, the finite element program is compiled with fully implicit integration algorithm which is return mapping algorithm of non-associative isotropic hardening Drucker-Prager criterion based on the elastic-plastic mechanics theory frame and the nonlinear finite element theory. Return mapping algorithm can avoid the drift phenomenon of the trial stress, and can achieve the accurate solution of the constitutive equation on the condition of the quasi-static deformation, a quadratic convergence rate when using the Newton-Raphson iteration scheme, higher accuracy and stability. To solve the problems of foundations in geotechnical engineering, displacements and stresses are calculated, and the evolution process of plastic zones is simulated. The ultimate bearing capacity of the analytical solution is compared with that from the numerical solution. The results demonstrate the superiority of the algorithm, and the correctness and the practicality of the program.