The poloidal harmonics of the toroidal normal modes of an unstable axisymmetric tokamak plasma are employed as basis functions for the minimisation of the 3D energy functional. This approach presents a natural extension of the perturbative method considered in [M.S. Anastopoulos Tzanis et al, Nuclear Fusion 59:126028, 2019]. This variational formulation is applied to the stability of tokamak plasmas subject to external non-axisymmetric magnetic fields. A comparison of the variational and perturbative methods shows that for D-shaped, high betaN plasmas, the coupling of normal modes becomes strong at experimentally relevant applied 3D fields, leading to violation of the assumptions that justify a perturbative analysis. The variational analysis employed here addresses strong coupling, minimising energy with respect to both toroidal and poloidal Fourier coefficients. In general, it is observed that ballooning unstable modes are further destabilised by the applied 3D fields and field-aligned localisation of the perturbation takes place, as local ballooning theory suggests. For D-shaped high betaN plasmas, relevant to experimental cases, it is observed that the existence of intermediate n unstable peeling-ballooning modes, where a maximum in the growth rate spectrum typically occurs, leads to a destabilising synergistic coupling that strongly degrades the stability of the 3D system.