This manuscript proposes a novel, efficient finite element solution technique for the computational simulation of cardiac electrophysiology. We apply a two-parameter model that is characterized through a fast action potential and a slow recovery variable. The former is introduced globally as a nodal degree of freedom, whereas the latter is treated locally as internal variable on the integration point level. This particular discretization is extremely efficient and highly modular since different cardiac cell models can be incorporated straightforwardly through only minor local modifications on the integration point level. In this manuscript, we illustrate the algorithm in terms of the Aliev-Panfilov model for cardiomyocytes. To ensure unconditional stability, a backward Euler scheme is applied to discretize the evolution equation for both the action potential and the recovery variable in time. To increase robustness and guarantee optimal quadratic convergence, we suggest an incremental iterative Newton-Raphson scheme and illustrate the consistent linearization of the weak form of the excitation problem. The proposed algorithm is illustrated by means of two- and three-dimensional examples of re-entrant spiral and scroll waves characteristic of cardiac arrhythmias in atrial and ventricular fibrillation.