This paper aims to investigate the differences in factor of safety (FS) and failure mechanism (FM) for spatially variable undrained soil slope between using finite element method (FEM) , finite difference method (FDM), and limit equilibrium method (LEM). The undrained shear strength of cohesive soil slope is modeled by a one-dimensional random field in the vertical direction. The FS and FM for a specific realization of random field are determined by SRT embedded in FEM- and FDM-based software (e.g., Phase2 6.0 and FLAC) and LEM, respectively. The comparative study has demonstrated that the bishop method (with circular failure surface) exhibits performance as fairly good as that of SRT both in FS and FM for the undrained slope cases where no preferable controlling surfaces such as hydraulic tension crack and inclined weak seams dominate the failure mechanism. It is, however, worthwhile to point out that unconservative FM is provided by the Bishop method from the aspect of failure consequence (i.e., the failure consequence indicated by the FM from the Bishop method is smaller than that from SRT). The rigorous LEM (e.g., M-P and Spencer method with noncircular failure surface) is not recommended in the stability analysis of spatially variable soil slopes before the local minima and failure to converge issues are fully addressed. The SRT in combination with FEM and/or FDM provides a rigorous and powerful tool and is highly preferable for slope reliability of spatially variable undrained slope.