The modeling of the controlled-source electromagnetic (CSEM) and single-well and crosswell electromagnetic (EM) configurations requires fine gridding to take into account the 3D nature of the geometries encountered in these applications that include geological structures with complicated shapes and exhibiting large variations in conductivities such as the seafloor bathymetry, the land topography, and targets with complex geometries and large contrasts in conductivities. Such problems significantly increase the computational cost of the conventional finite-difference (FD) approaches mainly due to the large condition numbers of the corresponding linear systems. To handle these problems, we employ a volume integral equation (IE) approach to arrive at an effective preconditioning operator for our FD solver. We refer to this new hybrid algorithm as the finite-difference integral equation method (FDIE). This FDIE preconditioning operator is divergence free and is based on a magnetic field formulation. Similar to the Lippman-Schwinger IE method, this scheme allows us to use a background elimination approach to reduce the computational domain, resulting in a smaller size stiffness matrix. Furthermore, it yields a linear system whose condition number is close to that of the conventional Lippman-Schwinger IE approach, significantly reducing the condition number of the stiffness matrix of the FD solver. Moreover, the FD framework allows us to substitute convolution operations by the inversion of banded matrices, which significantly reduces the computational cost per iteration of the hybrid method compared to the standard IE approaches. Also, well-established FD homogenization and optimal gridding algorithms make the FDIE more appropriate for the discretization of strongly inhomogeneous media. Some numerical studies are presented to illustrate the accuracy and effectiveness of the presented solver for CSEM, single-well, and crosswell EM applications.