A three-dimensional (3D) volume integral equation solution was refined and adapted to magnetotelluric (MT) modeling. The refinement, incorporating an integro-difference scheme, increases the accuracy somewhat without increasing the computer time. Utilizing the two symmetry planes for a plane wave source decreases the computer storage by a factor of 8 and greatly reduces the computer time. Convergence checks and comparisons with other solutions show that our results are valid. Because of space charges at resistivity boundaries, low-frequency 3D responses are much different from 1D and 2D responses. Hence 3D models are required for interpreting MT data in the complex geothermal environment