<p>This study investigates coupled effects of strength anisotropy and strain softening in natural saturated clays on undrained slope stability, together with the associated numerical challenges. Traditional finite element methods suffer mesh dependency from strain softening, requiring regularization techniques like the Cosserat continuum. However, for large-scale clay analyses, a systematic method to determine the internal length scale parameter (<i>l</i>_c) is lacking. Moreover, the lack of comparative studies and guidance on parameter selection for different anisotropic strength methods limits engineering applications. To address these issues, this study systematically compares three methods: the Casagrande approach, the microstructure tensor-based Pietruszczak–Mroz (PM) method, and a proposed Modified PM method. Their respective recommended test types are conventional triaxial, uniaxial, and plane strain compression. A unified anisotropic evolution equation is formulated by combining these methods with an exponential strain softening law and implemented in a Tresca-inscribed Mises model within the Cosserat continuum framework. A practical approach to determine <i>l</i>_c in large-scale clay conditions is also proposed. Numerical results show that when the deposition angle <i>α</i> = 0°, the normalized stability number <i>N</i>_s ranks as Casagrande &gt; Modified PM &gt; PM for <i>k</i> &lt; 1.0, and reverses for<i> k</i> &gt; 1.0. For typical cases (<i>α</i> = 0° and k ≤ 1.0), <i>N</i>_s decreases from isotropy to inherent anisotropy, evolving anisotropy, and inherent anisotropy with strain softening, with greater differences at gentler slopes. In the San Francisco Bay case, the proposed method, informed by test data and incorporating anisotropy and strain softening, yielded <i>F</i>_s = 1.06–1.07, which is closer to the actual failure condition (<i>F</i>_s ≤ 1.0) than results reported in previous studies, confirming the model’s engineering applicability.</p>