Abstract. We study a numerical method to compute probability density functions of solutions of stochastic differential equations. The method is sometimes called the numerical path integration method and has been shown to be fast and accurate in application oriented fields. In this paper we provide a rigorous analysis of the method that covers systems of equations with unbounded coefficients. Working in a natural space for densities, L 1 , we obtain stability, consistency, and new convergence results for the method, new well-posedness and semigroup generation results for the related Fokker-Planck-Kolmogorov equation, and a new and rigorous connection to the corresponding probability density functions for both the approximate and the exact problems. To prove the results we combine semigroup and PDE arguments in a new way that should be of independent interest.