The　meshless　numerical　manifold　method　(MNMM)　inherits　the　advantages　of　handling　continuous　and　discontinuous　problems　uniformly,　improving　interpolation　precision,　and　avoiding　linear　dependence　issues.　This　paper　refines　the　MNMM　method　for　slope　stability　analysis　to　overcome　the　dependence　and　sensitivity　of　the　strength　reduction　method　of　the　finite　element　method　on　the　number　of　finite　element　meshes　and　the　algorithm　for　solving　the　nonlinear　equations　when　calculating　slopes　with　soft　and　weak　interlayers.　First,　the　strength　reduction　method　(SRM)　was　incorporated　into　MNMM　to　create　the　numerical　model　called　MNMM-SRM.　Based　on　the　MNMM-SRM,　a　complementary　Mohr-Coulomb　(M-C)　control　equation　compliant　with　Koiter＇s　rule　of　multi-yielding　surfaces　was　established.　In　the　principal　stress　space,　the　subspace　tracking　method　was　used　to　track　the　eigenvalues　and　eigenvectors　of　the　stress　tensor,　thereby　addressing　the　corner　problem　of　M-C　multi-yield　surfaces.　Combining　the　k-means　clustering　algorithm　with　wavelet　transformation　and　employing　the　slope＇s　total　nodal　displacement　field　as　an　input　variable,　an　intelligent　and　efficient　critical　slip　surface　(CSS)　automatic　extraction　method　was　developed.　Three　typical　cases　were　used　to　validate　the　accuracy　of　the　numerical　model.　This　numerical　model　has　a　wide　range　of　potential　applications　in　slope　stability　research.