A new method to compute the dynamic response of periodic waveguides with localised nonlinearities is introduced and used to investigate the nonlinear shift of a band-edge mode in the bandgap of a locally resonant phononic structure. This nonlinear extension of the Wave Finite Element Method (WFEM) uses a finite-element discretisation of arbitrarily complex unit-cells, and leverages Floquet-Bloch theory to reduce the analysis of the entire waveguide to a state-vector of Bloch waves' amplitude. Higher harmonics generated by nonlinear effects are addressed using the Harmonic Balance Method and the nonlinear forces are evaluated via an alternating frequency-time procedure. The periodic response of the system is computed through a continuation scheme, taking the Bloch waves' amplitude as unknowns. The accuracy of the nonlinear WFEM is validated against standard FEM with Craig-Bampton reduction, demonstrating an 83% speedup in resolution time. Applying the method to a locally resonant metamaterial demonstrates that nonlinear effects can shift resonances from outside to inside bandgaps, resulting in high-amplitude, spatially localised vibrations where small amplitudes are expected from linear theory. The versatility and computational efficiency of this nonlinear dynamic simulation method should facilitate the study of complex metamaterials and civil engineering structures coupled with nonlinear interfaces or singularities.