We have developed a new 3D magnetotelluric modeling scheme in a mixed space-wavenumber domain. The modeling scheme is based on using a 2D Fourier transform along two horizontal directions to solve a vector-scalar potential formula derived from Maxwell’s equations based on the primary-secondary potential separation. The derived 1D governing equations in a mixed space-wavenumber domain are solved by using the finite-element method (FEM) together with a chasing method, and then the 2D inverse Fourier transform is used to recover the final solution of the electromagnetic (EM) fields in the 3D spatial domain. An iterative scheme is applied to approximate the true solution by repeating the previous steps because the governing equations cannot be solved directly due to an unusual primary-secondary potential field separation used. Nevertheless, the new method is capable of reducing the memory requirement and computational time in the mixed domain, and the 1D governing equations are highly parallel among different wavenumbers. For each of the 1D equations, the two- or four-node Gaussian quadrature rule can be used in both horizontal directions for Gauss fast Fourier transform. It is worth mentioning that the linear matrix equation to be solved is a fixed bandwidth system, and the chasing method is more efficient and convenient than solvers with preconditioners for the 1D matrix equations. The reliability and efficiency of the newly proposed method are verified with three synthetic 3D models by comparisons with a classical integral equation solution, an adaptive FEM solution, and a nonadaptive FEM solution. The proposed algorithm will be used in electrical resistivity tomography and controlled-source EM methods in future studies.