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 anisotropic diffusion equation in the presence of advection and time-independent as well as 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 our subsequent works, 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.