AbstractThis paper extends Eshelby’s problem of one inhomogeneity embedded in a homogeneous infinite domain to a bimaterial infinite domain. The equivalent inclusion method (EIM) was used to simulate the inhomogeneity of an inclusion with a polynomial eigenstrain. The fundamental solution of a point force in a bimaterial was used to formulate the domain integral over the inclusion. For a finite bimaterial domain, the boundary integral equation (BIE) takes into account the boundary responses by a single domain instead of utilizing the conventional multiregion BIE scheme. The EIM can be used similarly, and the elastic field can be obtained with tailorable accuracy based on the order of the polynomial eigenstrain. The algorithm is particularly suitable to simulate a defect in thin film–substrate systems or other similar bilayered materials. The stress concentration of a microvoid embedded in a bilayered solar panel was investigated. The size and location of the void relative to the interface exhibited considerable effects on the stress concentration factor. Numerical case studies demonstrated the effectiveness and accuracy of the algorithm, and parametric studies showed the boundary effects on the stress concentration of a microvoid in a finite bimaterial under a uniform far-field strain.