This work introduces a novel, unconditionally stable and fully coupled finite element method for the bidomain system of equations of cardiac electrophysiology. The transmembrane potential phi(i) - phi(e) and the extracellular potential phi(e) are treated as independent variables. To this end, the respective reaction-diffusion equations are recast into weak forms via a conventional isoparametric Galerkin approach. The resultant nonlinear set of residual equations is consistently linearised. The method results in a symmetric set of equations, which reduces the computational time significantly compared to the conventional solution algorithms. The proposed method is inherently modular and can be combined with phenomenological or ionic models across the cell membrane. The efficiency of the method and the comparison of its computational cost with respect to the simplified monodomain models are demonstrated through representative numerical examples.