  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.
