Abstract: In this paper, a method dividing a big element into smaller meshes is adopted to judge the ranges of plastic zones. The analytical function series is unified in both elastic and plastic zones. The method of modified stiffness increment is used in iteration,and the plastic coefficient matrix and plastic stiffness matrix are given in derivation. The element stresses in the mesh are regarded as constants that makes the solution process simpler.The program is developed to analyze the stability of three-dimentional underground opening, and the stability of circular and elliptic openings are also analyzed.