The bidomain and monodomain models are among the most widely used mathematical models to describe cardiac electrophysiology. They take the form of multi-scale reaction-diffusion partial differential equations that couple the dynamic behaviour on the cellular scale with that on the tissue scale. The systems of differential equations associated with these models are large and strongly non-linear, but they also have a distinct structure due to their multi-scale nature. For these reasons, numerical solutions to these systems are often found via operator-splitting methods. The focus of this presentation is on operator-splitting methods with order higher than two that, according to the Sheng--Suzuki theorem, require backward time integration and historically have been considered unstable for solving deterministic parabolic systems. The stability and performance of operator-splitting methods of up to order four to solve the bidomain and monodomain models are demonstrated on several examples arising in the field of cardiovascular modeling.