We developed an algorithm to simulate magnetotelluric responses for three-dimensional general electrical anisotropic structures with finite element (FE) method. From the Maxwell's equations, the curl-curl equation in terms of the electric vector potential and the magnetic scalar potential is obtained by decomposing the magnetic field into the incident field and the resultant field. For this method, the only unknown variable in the air is the magnetic scalar potential. Galerkin weighted residual method is adopted to generate the variational equations. Then we made a comparison with the edge-based FE method which is in terms of electric fields to validate the accuracy of this method. Later, two types of anisotropic models are studied and the results are analyzed in detail. Finally, based on these analyses, two main meaningful conclusions are obtained.