This study first reviews the numerical manifold method (NMM) which possesses some advantages over the traditional limit equilibrium methods (LEMs) in calculating the factors of safety (Fs) of the slopes. Then, with regard to a trial slip surface (TSS), associated stress fields reproduced by NMM as well as the enhanced limit equilibrium method are combined to compute Fs. In order to search for the potential critical slip surface (CSS), the MAX-MIN ant colony optimization algorithm (MMACOA), one of the best performing algorithms for some optimization problems, is adopted. Procedures to obtain Fs in conjunction with the potential CSS are described. Finally, the proposed numerical model and traditional methods are compared with stability analysis of three typical slopes. The numerical results show that Fs and CSSs of the slopes can be accurately calculated with the proposed model.