We present a finite volume (FV) algorithm for the magnetotelluric (MT) forward modeling in 3-D 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 employed 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 1-D anisotropic model. For this model, the FV code provides accurate results, and the computational cost is reasonable. The iterative solvers QMR and BiCGStab preconditioned with the electromagnetic potential (A-ø) system exhibit 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 seconds with multiple threads involved in the computation, and the memory usage is only 11.6 GB in the “Out-of-Core” mode. We further calculate MT responses of a 3-D model with dipping and horizontal anisotropy, respectively. The results suggest that electrical anisotropy can have significant influence on the MT response.