A numerical technique has been developed to solve the three-dimensional (3-D) potential distribution about a point source of current located in or on the surface of a half-space containing an arbitrary 3-D conductivity distribution. Self-adjoint difference equations are obtained for Poisson's equation using finite-difference approximations in conjunction with an elemental volume discretization of the lower half-space. Potential distribution at all points in the set defining the subsurface are simultaneously solved for multiple point sources of current. Accurate and stable solutions are obtained using full, banded, Cholesky decomposition of the capacitance matrix as well as the recently developed incomplete Cholesky-conjugate gradient iterative method. A comparison of the 2-D and 3-D simple block-shaped models, for the collinear dipole-dipole array, indicates substantially lower anomaly indices for inhomogeneities of finite strike-extent. In general, the strike-extents of inhomogeneities have to be approximately 10 times the dipole lengths before the response becomes 2-D. The saturation effect with increasing conductivity contrasts appears sooner for the 3-D conductive inhomogeneities than for corresponding models with infinite strike-lengths.