Recent proceedings in radiation belt studies have proposed new requirements for numerical methods to solve the kinetic equations involved. In this article, we present a numerical solver that can solve the general form of the radiation belt Fokker–Planck equation and Boltzmann equation in arbitrarily provided coordinate systems and with user-specified boundary geometry, boundary conditions, and equation terms. The solver is based upon the mathematical theory of stochastic differential equations, whose computational accuracy and efficiency are greatly enhanced by specially designed adaptive algorithms and a variance reduction technique. The versatility and robustness of the solver are exhibited in four example problems. The solver applies to a wide spectrum of radiation belt modeling problems, including the ones featuring non-diffusive particle transport such as that arising from nonlinear wave–particle interactions.