Abstract: Equivalent inclusion method is a convenient tool for modeling the stress field of materials embedded with inhomogeneities. However, its analytical applications are limited to elliptical shaped inhomogeneity cases. In this work, a new powerful and versatile numerical extension to the original equivalent inclusion method, called numerical equivalent inclusion method, was proposed. The fundamental theory of the numerical equivalent inclusion method was introduced, and an implementation method, conjugate gradient method, was presented to solve the linear equation group of the equivalency condition of the new method. The method can be applied to 2D inhomogeneity problems with arbitrary shape through a handy numerical discretization. Benchmark comparisons with the analytical results for an elliptical inhomogeneity model illustrated the accuracy of the proposed solution method. The efficiency and convergence of the numerical equivalent inclusion method were also discussed in detail, the results show that the proposed method implemented by conjugate gradient method has significant advantages in efficiency compared with that by Gaussian elimination method, and can maintain the accuracy of the results as well. A half-elliptical inhomogeneity model and a kind of zirconia/alumina coextruded composites were utilized to demonstrate the capability of the new method on solving arbitrarily shaped inhomogeneity problems.