  We present a new numerical algorithm to solve the recently derived equations
of two-moment cosmic ray hydrodynamics (CRHD). The algorithm is implemented as
a module in the moving mesh Arepo code. Therein, the anisotropic transport of
cosmic rays (CRs) along magnetic field lines is discretised using a
path-conservative finite volume method on the unstructured time-dependent
Voronoi mesh of Arepo. The interaction of CRs and gyroresonant Alfv\'en waves
is described by short-timescale source terms in the CRHD equations. We employ a
custom-made semi-implicit adaptive time stepping source term integrator to
accurately integrate this interaction on the small light-crossing time of the
anisotropic transport step. Both the transport and the source term integration
step are separated from the evolution of the magneto-hydrodynamical equations
using an operator split approach. The new algorithm is tested with a variety of
test problems, including shock tubes, a perpendicular magnetised discontinuity,
the hydrodynamic response to a CR overpressure, CR acceleration of a warm
cloud, and a CR blast wave, which demonstrate that the coupling between CR and
magneto-hydrodynamics is robust and accurate. We demonstrate the numerical
convergence of the presented scheme using new linear and non-linear analytic
solutions.
