The cardiac Purkinje fibres provide an important stimulus to the coordinated contraction of the heart. We present a numerical algorithm for the solution of electrophysiology problems on the Purkinje network that is efficient enough to be used on realistic networks with physiologically detailed membrane models. The algorithm is based on operator splitting and is provided with three different implementations: pure CPU, hybrid CPU/GPU, and pure GPU. Compared to our previous work based on the model of Vigmond et al., we modify the explicit gap junction term at network bifurcations in order to improve its mathematical consistency. Due to this improved consistency of the model, we are able to perform a convergence study against analytical solutions and verify that all three implementations produce equivalent convergence rates. Finally, we compare the efficiency of all three implementations on Purkinje networks of increasing spatial resolution using membrane models of increasing complexity. Both hybrid and pure-GPU implementations outperform the pure-CPU implementation, but their relative performance difference depends on the size of the Purkinje network and the complexity of the membrane model used. Copyright c 2015 John Wiley & Sons, Ltd.