Strain softening problem in geotechnical engineering and the difficult solution problem of the finite element numerical calculation due to the negative tangent stiffness of strain softening model are studied. An elastoplastic strain softening constitutive model of rock is established based on the Drucker-Prager strength criteria. A fully implicit return mapping algorithm which has characteristics of the unconditional stability and precision is used to solve the constitutive equation, and how the programmed constitutive model to be solved is discussed in detail. Then, the shortcomings of low efficiency of the arc-length method in judging stiffness matrix is considered, Newton-Raphson scheme and arc-length method (NR-AL method) are combined to iteratively solve the calculation of elastoplastic incremental finite element equations. Namely Newton-Raphson scheme is used before the structure reaching the limit load, and when the structure is close to the limit load, turning to the arc-length method, so that the structure can go over the peak point into the softening phase until destruction. NR-AL method has the advantages in the iterative solution. A program of the built strain softening model and elastoplastic incremental finite element to solve the constitutive equation for the iterative process is compiled using C++ language. The program is applied to numerical calculation, and the stress-strain curves of the idealized elastoplastic model, strain softening and strain hardening model based on the Drucker-Prager strength criteria are comparatively analyzed. The results show that the strain softening constitutive model can simulate the characteristics the post-peak softening of rock material well, and it can reveal the features of the post-peak strain softening and failure mechanism. NR-AL method can solve the negative stiffness problem caused by strain softening and also overcome the shortcomings of low efficiency in judging stiffness matrix using the arc-length only.