Numerical Limit Analysis (NLA), although little known and used, is a powerful technique to efficiently determine the stability conditions in slopes. In this paper, we extend the application of the NLA by including a probabilistic framework to capture the spatial variability of soil shear strength parameters. To achieve this, stationary, unconditional, multivariate, cross-correlated random fields were implemented on an existing, in-house developed computer code, employing both conventional and stepwise Covariance Decomposition Methods (CMD). In order to verify the implemented code in terms of efficiency and robustness, 2D and 3D problems were analyzed. The paper demonstrates that NLA is a robust alternative to standard methods such as the limit equilibrium and finite element method with the strength reduction technique. Emphasis is placed on the significance of incorporating autocorrelation and cross-correlation to realistically represent spatial variation in parameters and address uncertainties associated with soil shear strength.