  We present the implementation of an implicit-explicit (IMEX) Runge-Kutta
numerical scheme for general relativistic hydrodynamics coupled to an optically
thick radiation field in two existing GR-hydrodynamics codes. We argue that the
necessity of such an improvement arises naturally in astrophysically relevant
regimes where the optical thickness is high as the equations become stiff. By
performing several 1D tests we verify the codes' new ability to deal with this
stiffness and show consistency. Then, still in 1D, we compute a luminosity
versus accretion rate diagram for the setup of spherical accretion onto a
Schwarzschild black hole and find good agreement with previous work. Lastly, we
revisit the supersonic Bondi Hoyle Lyttleton (BHL) accretion in 2D where we can
now present simulations of realistic temperatures, down to T~10^6 K. Here we
find that radiation pressure plays an important role, but also that these
highly dynamical set-ups push our approximate treatment towards the limit of
physical applicability. The main features of radiation hydrodynamics BHL flows
manifest as (i) an effective adiabatic index approaching gamma_effective ~ 4/3;
(ii) accretion rates two orders of magnitude lower than without radiation
pressure; (iii) luminosity estimates around the Eddington limit, hence with an
overall radiative efficiency as small as eta ~ 10^{-2}; (iv) strong departures
from thermal equilibrium in shocked regions; (v) no appearance of the flip-flop
instability. We conclude that the current optically thick approximation to the
radiation transfer does give physically substantial improvements over the pure
hydro also in set-ups departing from equilibrium, and, once accompanied by an
optically thin treatment, is likely to provide a fundamental tool for
investigating accretion flows in a large variety of astrophysical systems.
