This study performed a steady-state thermal analysis of ellipsoidal inhomogeneities embedded in a bimaterial or semi-infinite domain. Because of the Green’s function of the bimaterial, the continuity conditions for the temperature and heat flux of the interface were mathematically involved. Through proper modification of the material properties, the Green’s function of the bimaterial could be reduced to semi-infinite and infinite. This study proposed the use of Eshelby’s equivalent inclusion method to simulate ellipsoidal inhomogeneities, which replaced the inhomogeneity of the matrix material containing a continuously distributed polynomial-form eigen-temperature gradient field. Based on the analytical ellipsoidal integral with polynomial density functions, the disturbance caused by inhomogeneities was analytically evaluated using the domain integrals of the eigentemperature gradient and Green’s function of the bimaterial. The eigen-temperature gradient for each inhomogeneity was described by a Taylor series expanded at the geometric center, which could then be evaluated by solving the equivalent heat flux conditions. A fi- nite element method (FEM) was used to verify the accuracy of the semi-analytical method.The study showed that mesh-free solutions to multiple ellipsoidal inhomogeneities in the bimaterial/semi-infinite domain could be achieved.