Material nonlinearity and contact nonlinearity are prevalent in engineering structure service, and their induced mechanical response deviations, structural failure risks and coupling amplification effects seriously threaten service safety and lifespan. To address the dual nonlinearity of heterogeneous materials and contact, this study developed a numerical solution algorithm for singular integral equations of half-plane contact problems with impurities based on the Numerical Equivalent Inclusion Method (NEIM). The algorithm decomposes the heterogeneous material contact problem into an equivalent inclusion problem with unknown eigenstrains and a homogeneous material contact problem. Compared with the finite element method, it only requires mesh discretization of the contact zone and impurity region, significantly reducing computational memory and preprocessing burden. Verified by finite element results, the solution exhibits high accuracy; systematic analysis reveals the influence of elastic modulus ratios between impurities and matrix on contact stress distribution and internal stress fields. This work provides an efficient numerical method for engineering strength problems involving material-contact dual nonlinearity, such as gear contact fatigue.