Recent proceedings in the 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 radiation belt Fokker-Planck equation and Boltzmann equation in arbitrarily provided coordinate systems, and with user-specified boundary geometry and boundary conditions. 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 variance reduction technique. The versatility and robustness of the solver is exhibited in three example problems. The solver applies to a wide spectrum of radiation belt modeling problems, including the ones featuring nonlinear wave-particle interactions.