Heat transfer in magnetically confined plasmas is characterized by extremely high anisotropic diffusion phenomena. At the core of a magnetized plasma, the heat conductivity coefficients in the parallel and perpendicular directions of the induction field can be very different. Their ratio can exceed 1.e8 and the pollution by purely numerical errors can make very diffcult the simulation of the heat transport in the perpendicular direction. Standard numerical methods, generally used in the discretization of classical diffusion problems, are rather inefficient. The present paper analyses a finite element approach for the solution of a highly anisotropic diffusion equation. Two families of finite elements of class C1, namely bi-cubic Hemite-Bézier and reduced cubic Hsieh-Clough-Tocher finite elements, are compared. Their performances are tested numerically, for various ratios of the diffusion coefficients, on different mesh configurations, even aligned with the induction field. The time stepping is realised by an implicit high-order Gear finite difference scheme. An example of reduced model is also provided in order to comment on some obtained results.