A new code for solving radiation belt diffusion equations has been developed and applied to the 2-D bounce-averaged energy pitch angle quasi-linear diffusion equation. The code uses Monte Carlo methods to solve Itô stochastic differential equations (SDEs) which are mathematically equivalent to radiation belt diffusion equations. We show that our SDE code solves the diffusion equation with off-diagonal diffusion coefficients in contrast to standard finite difference codes which are generally unstable when off-diagonal diffusion coefficients are included. Our results are in excellent agreement with previous results. We have also investigated effects of assuming purely parallel propagating electromagnetic waves when calculating the diffusion coefficients and find that this assumption leads to errors of more than an order of magnitude in flux at some equatorial pitch angles for the specific chorus wave model we use. Further work is needed to investigate the sensitivity of our results to the wave model parameters. Generalization of the method to 3-D is straightforward, thus making this method a very promising new way to investigate the relative roles of pitch angle, energy, and radial diffusion in radiation belt dynamics.