  A numerical method for integrating the equations describing a dynamically
coupled system made of a fluid and cosmic-rays is developed. In smooth flows
the effect of CR pressure is accounted for by modification of the
characteristic equations and the energy exchange between cosmic-rays and the
fluid, due to diffusive processes in configuration and momentum space, is
modeled with a flux conserving method. Provided the shock acceleration
efficiency as a function of the upstream conditions and shock Mach number, we
show that the Riemann solver can be modified to take into account the
cosmic-ray mediation without having to resolve the cosmic-ray induced
substructure. Shocks are advanced with Glimm's method which preserves their
discontinuous character without any smearing, thus allowing to maintain
self-consistency in the shock solutions. In smooth flows either Glimm's or a
higher order Godunov's method can be applied, with the latter producing better
results when approximations are introduced in the Riemann solver.
