Slope stability analysis is a geotechnical engineering problem characterized by many sources of uncertainty. Some of them are connected to the uncertainties of soil properties involved in the analysis. In this paper, a numerical procedure for integrating commercial finite difference method into probabilistic analysis of slope stability is presented. Since the limit-state function cannot be expressed in an explicit form, the ANN-based response surface method is adopted to approximate the limit-state function and the first-, second-order reliability method and the Monte Carlo simulation technique are used to calculate the probability of failure. Probabilistic stability assessments for a hypothetical two-layer slope and the Sugar Creek embankment were performed to verify the application potential to the slope stability problems. The examples show the successful implementation and the possibility of the extension of the proposed procedure to the variety of geotechnical engineering problems.