Abstract: Cosmic ray (CR) protons are an important component in many astrophysical systems. Processes like CR injection, cooling, adiabatic changes as well as active CR transport through the medium strongly modify the CR momentum distribution and have to be taken into account in hydrodynamical simulations. We present an efficient novel numerical scheme to accurately compute the evolution of the particle distribution function by solving the Fokker-Planck equation with a low number of spectral bins (10 - 20), which is required to include a full spectrum for every computational fluid element. The distribution function is represented by piecewise power laws and is not forced to be continuous, which enables an optimal representation of the spectrum. The Fokker-Planck equation is solved with a two-moment approach evolving the CR number and energy density. The low numerical diffusion of the scheme reduces the numerical errors by orders of magnitude in comparison to classical schemes with piecewise constant spectral representations. With this method not only the spectral evolution of CRs can be computed accurately in magnetohydrodynamic simulations but also their dynamical impact as well as CR ionisation. This allows for more accurate models for astrophysical plasmas, like the interstellar medium, and direct comparisons with observations.