An algorithm is described to efficiently compute the cumulative distribution and probability density functions of the diffusion process (Ratcliff, 1978) with trial-to-trial variability in mean drift rate, starting point, and residual reaction time. Some, but not all, of the integrals appearing in the model’s equations have closed-form solutions, and thus we can avoid computationally expensive numerical approximations. Depending on the number of quadrature nodes used for the remaining numerical integrations, the final algorithm is at least 10 times faster than a classical algorithm using only numerical integration, and the accuracy is slightly higher. Next, we discuss some special cases with an alternative distribution for the residual reaction time or with fewer than three parameters exhibiting trialto-trial variability.