It is important to study the responses of the magnetotelluric (MT) method in anisotropic media. However, MT anisotropy research has focused mainly on one-dimensional (1D) and two-dimensional (2D) solutions. Therefore, we developed a three-dimensional (3D) finite element (FE) algorithm for MT modelling in anisotropic media. This approach is based on the weak formulation of the governing Maxwell equations using Coulomb-gauged potentials. The node-based FE method is adopted here, and the values of the coefficient matrixes are obtained with hexahedral meshes. To validate the correctness and accuracy of this method, its results are compared with previous solutions for a 2D anisotropic model and a 3D arbitrary anisotropic model, respectively. Different solvers with different preconditioners are tested, and the results show that the quasi-minimum residual method with the incomplete LU preconditioner is more stable and faster compared with the other schemes. We then studied a 3D anisotropic model in three different conditions, and analysed the results in detail. Finally, three main conclusions are obtained: the xy- and yy-mode apparent resistivities remain almost unchanged if a principal conductivity is in the x-direction; the yx- and xx-mode apparent resistivities remain almost the same if a principal conductivity is in the y-direction; a principal conductivity in the z-direction has almost no influence on apparent resistivities.