The determination of three-dimensional (3D) critical slip surface and its associated factor of safety (FoS) of slopes are essential to their stability analysis. However, it is difficult to obtain 3D slip surface conveniently and reliably by using existing analytical methods. In this paper, the FoS of slopes is calculated based on the stresses from the finite element method (FEM). A new method is developed to search 3D critical slip surface of slope by constructing an extended ellipsoidal slip surface and adopting self-adaptive differential evolution (SaDE) algorithm in combination with the opposite-movement scheme. The extended ellipsoidal slip surface is easy to construct and flexible to search, and meanwhile it is more coincident with the actuals than spherical and ellipsoidal slip surface. In addition, the whole calculating process is explicit without the convergence problem. The application of the developed method is verified by re-analyzing benchmark slope stability examples from the literature. The results obtained in examples prove the SaDE algorithm has high efficiency and global optimization ability, which are also in good agreement with limit equilibrium and strength reduction methods. It is also demonstrated that this method is suitable to be applied in stability analysis of slopes under complex lithology and environmental conditions by analyzing actual engineering slopes.