We present Sapphire++, an open-source code designed to numerically solve the Vlasov-Fokker-Planck equation for astrophysical applications. Sapphire++ employs a numerical algorithm based on a spherical harmonic expansion of the distribution function, expressing the Vlasov-Fokker-Planck equation as a system of partial differential equations governing the evolution of the expansion coefficients. The code utilises the discontinuous Galerkin method in conjunction with implicit and explicit time stepping methods to compute these coefficients, providing significant flexibility in its choice of spatial and temporal accuracy. We showcase the code's validity using examples. In particular, we simulate the acceleration of test particles at a parallel shock and compare the results to analytical predictions. The Sapphire++ code (https://github.com/sapphirepp/sapphirepp) is available as a free and open-source tool for the community.

25 pages, 11 figures, Published in Journal of Computational Physics