A finite element method based on bicubic Bézier surfaces is applied to the simulation of MHD instabilities relevant to magnetically confined fusion. The major advantage of the new technique is that it allows a natural way to implement mesh refinement strategy, which is not supported by a pure Hermite formulation. Compared to a Lagrangian formulation the number of degrees of freedom is significantly reduced. The use of an isoparametric representation of the space coordinates allows an accurate alignment of the finite elements to the magnetic field line geometry in a tokamak plasma. The Bézier finite elements have been implemented in a MHD code using the non-linear reduced MHD model in toroidal geometry. As an illustration, results for Soloviev equilibrium and time-dependent current-hole computations are presented and discussed.