We have developed a finite volume (FV) algorithm for magnetotelluric (MT) forward modeling in 3D conductivity structures with general anisotropy. The electric and magnetic fields are discretized on a conventional staggered grid, which cannot directly address the full-tensor conductivity. To overcome this difficulty, an interpolation scheme is used to average different components of the electric field to the same position. We formulate the algorithm in pure matrix form and implement it in a new language, Julia, making the programming process highly efficient and leading to a code with excellent readability, maintainability, and extendability. The validity of the FV Julia code is demonstrated using a layered 1D anisotropic model. For this model, the FV code provides accurate results, and the computational cost is reasonable. Being preconditioned with the electromagnetic potential () system, the iterative solvers including quasi-minimal residual and biconjugate gradient stabilized exhibit a good convergence rate for a wide range of periods. The direct solvers MUMPS and PARDISO are highly efficient for small model sizes. For a relatively large model size with 2.18 millions unknowns, the linear system of one period can be solved by MUMPS within 360 s with multiple threads involved in the computation, and the memory usage is only 11.6 GB in the “out-of-core” mode. We further calculated MT responses of a 3D model with dipping and horizontal anisotropy, respectively. The results suggest that the electrical anisotropy can have significant influence on the MT response.