This paper develops a numerical method for solving multiple three-dimensional inhomogeneous inclusions of arbitrary shape in an isotropic half space under external loading. The method considers interactions between all the inhomogeneous inclusions and thus could provide an accurate stress field for the analysis of material strength and reliability. In the method, the inhomogeneous inclusions are first broken up into small cuboidal elements, which each are then treated as cuboidal homogeneous inclusions with initial eigenstrains plus unknown equivalent eigenstrains using Eshelby’s equivalent inclusion method. The unknown equivalent eigenstrains are introduced to represent the material dissimilarity of the inhomogeneous inclusions, their interactions and their response to external loading, and determined by solving a set of simultaneous constitutive equations established for each equivalent cuboidal inclusion. The method is validated by the finite element method and then applied to investigate a cavity-contained inhomogeneous inclusion and a stringer/cluster of inhomogeneities near a half-space surface. This solution may have potentially significant application in addressing challenging material science and engineering problems concerning inelastic deformation and material dissimilarity.