Methods for solving elastic wave scattering problems in three dimensions (3D) with multiple inhomogeneities are discussed. A problem of homogeneous, isotropic elastic defects in an otherwise homogeneous, isotropic elastic full‐space is formulated as a boundary integral equation. This equation is solved by discretizing the surface of each scatterer in a fashion known as the boundary element method. The resulting matrix equation may be solved in a fully implicit manner, but an implicit‐iterative method is more efficient. With this hybrid method, a portion of the nonsingular integral operator is expanded in a Neumann series. Terms in this series correspond physically to Nth‐order Born approximations of the scatterers’ interaction. The relative advantage of this hybrid scheme depends on the number of iterations required. Except for closely situated strong scatterers, terms beyond the first few orders are not significant and thus the method can be quite advantageous. When the separation is large, another approximate method which ignores the evanescent portion of the near scattered field and further neglects the curvature of scattered waves is appropriate. Results from the converged, implicit‐iterative approach are compared with this far‐field approximation for many situations involving spheroidal voids. The validity of this approximation is explored in the near field.