We present a new, accurate, high-performance, wide-band three-dimensional (3-D) solver for the electromagnetic (EM) field scattering problem in an isotropic earth. The solver relates to those based on the volume integral equation (IE) approach and exploits a modified Neumann series (MNS) technique to solve Maxwell's equations. The solver allows for the conduction, polarization and displacement currents to be taken into account and admits for 3-D earth excitation by arbitrary electric or/and magnetic sources. We estimate the solver efficiency for scatterers discretized into N-x x N-y x N-z prisms, where it requires only about 6N(x)N(y)N(z)(log(2)(2N(x))log(2)(2N(y)) + 6N(z)) multiplications to get one term of the MNS expansion and about 200 NxNyNz2 bytes of memory. Our experience show that the number of terms N which are to be summed up to get the solution to 1% accuracy doesn't exceed fifty for the models with the conductivity contrast of up to 100. We demonstrate the solver versatility for magnetotellurics (MT) and controlled-source simulations. EM fields arising from a 3-D model with two high-contrast thin layers residing in layered earth were simulated due to a 10 Hz electric dipole located at the surface. When the layers were discretized into 16,384 prisms our code on a Pentium-100 MHz took T similar to 58 minutes, M similar to 7 Mbytes and N similar to 280. We also modeled the 0.1 Hz and 0.01 Hz MT responses within 3-D model with 1 Ohm.m and 100 Ohm.m blocks. When the blocks were discretized into 8,000 prisms the code took T similar to 5 minutes, M similar to 8 Mbytes, and N similar to 25. Finally fields for a crosswell model including a 3-D conducting target were simulated for 0.1 kHz and 10 kHz electric and magnetic dipoles in the wellbores. While the target was discretized into 6,250 prisms the code took T similar to 16 minutes, M similar to 13 Mbytes, and N similar to 24. All simulations showed from very good to excellent agreement with those of the other 3-D solvers.