A finite difference algorithm for 3D modeling of magnetotelluric data is presented. The algorithm, MT_3D_EA, (Magnetotelluric 3-Dimensional Eigenmode Approach) uses eigenmodes for solving the relevant partial differential equation wherein only a small subset of eigenvalues and corresponding eigenvectors are required to achieve satisfactory results. The requisite small subset (pre-specified number) of eigenmodes are obtained by using shift and invert implementation of Implicitly Restarted Arnoldi Method (IRAM). It is analyzed that only 10-15% smallest eigenvalues and corresponding eigenvectors are sufficient to obtain the acceptable solution. Using this approach one can find multi-frequency solution almost at the computational cost of single frequency solution. 3D results for synthetic models taken from COMMEMI project report and simulation results for simplified geoelectrical model of Garhwal Himalaya based on field data are presented as experiment design.