As the international tokamak ITER is being built, non-linear MHD simulations are playing an essential role in active research, understanding, and prediction of tokamak plasmas for the realisation of a fusion power plant. The development of MHD codes like JOREK is a key aspect of this research effort, and provides invaluable insight into the plasma stability and the control of global and localised plasma events, like Edge-Localised-Mode and disruptions. In this paper, we present an operational implementation of a new, generalised formulation of Bezier finite-elements applied to the JOREK code, a significant advancement from the previously C1-continuous bi-cubic Bezier elements. This new mathematical method enables any polynomial order of Bezier elements, with a guarantee of C-continuity at the level of (n−1)/2, where n is the order of the Bezier polynomials. The generalised method is defined, and a rigorous mathematical proof is provided for the C-continuity requirement. Key details on the code implementation are mentioned, together with a suite of tests to demonstrate the mathematical reliability of the finite-element method, as well as the practical usability for typical non-linear tokamak MHD simulations. These tests show very clearly the many advantages of this new method. A demonstration for a state-of-the-art simulation of an ITER Edge-Localised-Mode, with complex grid geometry, finalises the study