  We present a new approach for stably evolving general relativistic
magnetohydrodynamic (GRMHD) simulations in regions where the magnetization
$\sigma=b^2/\rho c^2$ becomes large. GRMHD codes typically struggle to evolve
plasma above $\sigma\approx100$ in simulations of black hole accretion. To
ensure stability, GRMHD codes will inject mass density artificially to the
simulation as necessary to keep the magnetization below a ceiling value
$\sigma_{\rm max}$. We propose an alternative approach where the simulation
transitions to solving the equations of general relativistic force-free
electrodynamics (GRFFE) above a magnetization $\sigma_{\rm trans}$. We augment
the GRFFE equations in the highly magnetized region with approximate equations
to evolve the decoupled field-parallel velocity, plasma energy density, and
plasma mass density. Our hybrid scheme is explicit and easily added to the
framework of standard-volume GRMHD codes. We present a variety of tests of our
method, implemented in the GRMHD code KORAL, and we show first results from a
3D hybrid GRMHD+GRFFE simulation of a magnetically arrested disc (MAD) around a
spinning black hole. Our hybrid MAD simulation closely matches the average
properties of a standard GRMHD MAD simulation with the same initial conditions
in low magnetization regions, but it achieves a magnetization
$\sigma\approx10^6$ in the evacuated jet funnel. We present simulated
horizon-scale images of both simulations at 230 GHz with the black hole mass
and accretion rate matched to M87*. Images from the hybrid simulation are less
affected by the choice of magnetization cutoff $\sigma_{\rm cut}$ imposed in
radiative transfer than images from the standard GRMHD simulation.
