A refined numerical method for the inclusion problem involving known, but otherwise arbitrary shaped multiple regions with eigenstrains is advanced in this paper. The newly advanced algorithm, employing the mirror image method, derives its effectiveness from implementation of a hybrid multi-dimensional algorithm, based on convolution and correlation theorems. A good validation for various regularly shaped inclusions containing uniform eigenstrains is found. The method is subsequently applied to predict a critical interaction distance between two neighboring spherical inclusions, containing uniform or non-uniform dilatational eigenstrains.DOI: http://dx.doi.org/10.5755/j01.mech.19.3.4657