We develop and analyze a method, density tracking by quadrature (DTQ), to compute the probability density function of the solution of a stochastic differential equation. The derivation of the method begins with the discretization in time of the stochastic differential equation, resulting in a discrete-time Markov chain with continuous state space. At each time step, DTQ applies quadrature to solve the Chapman-Kolmogorov equation for this Markov chain. In this paper, we focus on a particular case of the DTQ method that arises from applying the Euler-Maruyama method in time and the trapezoidal quadrature rule in space. Our main result establishes that the density computed by DTQ converges in $L^1$ to both the exact density of the Markov chain (with exponential convergence rate), and to the exact density of the stochastic differential equation (with first-order convergence rate). We establish a Chernoff bound that implies convergence of a domain-truncated version of DTQ. We carry out numerical tests to show that the empirical performance of DTQ matches theoretical results, and also to demonstrate that DTQ can compute densities several times faster than a Fokker-Planck solver, for the same level of error.