Computational models of cardiac electrophysiology require the integration of large and complex systems of ordinary differential equations. This poses special challenges on a GPU architecture owing to the large numbers of state variables required to describe the detailed biophysics of the ion channels in the cardiac muscle cells. We describe a system for generating efficient CUDA C from high level mathematical expressions and their solution within large three-dimensional, anatomically realistic multi-scale models of whole heart electrical activation. Significant performance improvements over CPU clusters are presented.