The bidomain model is a popular model for simulating electrical activity in cardiac tissue. It is a continuum-based model consisting of non-linear ordinary differential equations (ODEs) describing spatially averaged cellular reactions and a system of partial differential equations (PDEs) describing electrodiffusion on tissue level. Because of this multi-scale, ODE/PDE structure of the model, operator-splitting methods that treat the ODEs and PDEs in separate steps are natural candidates as numerical solution methods. Second-order methods can generally be expected to be more effective than first-order methods under normal accuracy requirements. However, the simplest and the most commonly applied second-order method for the PDE step, the Crank-Nicolson (CN) method, may generate unphysical oscillations. In this paper, we investigate the performance of a two-stage, L-stable singly diagonally implicit Runge-Kutta method for solving the PDEs of the bidomain model. Numerical experiments show that the enhanced stability property of this method leads to more physically realistic numerical simulations compared to both the CN and backward Euler methods.