  Anisotropic diffusion is imperative in understanding cosmic ray diffusion
across the Galaxy, the heliosphere, and the interplay of cosmic rays with the
Galactic magnetic field. This diffusion term contributes to the highly stiff
nature of the cosmic ray transport equation. To conduct numerical simulations
of time-dependent cosmic ray transport, implicit integrators (namely,
Crank-Nicolson (CN)) have been traditionally favoured over the CFL-bound
explicit integrators in order to be able to take large step sizes. We propose
exponential methods to treat the linear anisotropc diffusion equation in the
presence of advection and time-independent and time-dependent sources. These
methods allow us to take even larger step sizes that can substantially speed-up
the simulations whilst generating highly accurate solutions. In or subsequent
work, we will use these exponential solvers in the Picard code to study
anisotropic cosmic ray diffusion and we will consider additional physical
processes such as continuous momentum losses and reacceleration.
