The recent development of semi-analytical methods (SAM) has led to numerous improvements in their capabilities in terms of phenomena that can be accounted for and numerical efficiency. They now allow to perform fast and robust simulations of contact between inelastic—with either elastic–plastic or viscoelastic behavior—and anisotropic or heterogeneous materials. All effects may be combined, with either coating, inclusions, cavities, or fibers as inhomogeneities. The coupling between local and global scales remains numerically difficult. A framework is proposed here for contact problems considering the effect of elastic heterogeneities within an elastic–plastic matrix. The mutual interactions among heterogeneities and their surrounding plastic zone as well as the interactions between them and the contact surface through which the load is transmitted should be accounted for. These couplings are outside the validity domain of the Eshelby’s equivalent inclusion method (EIM) that assumes a uniform stress field in an infinite space far from the inhomogeneity. In the presence of heterogeneities close to the surface or located at the Hertzian depth, the yield stress can be reached locally due to the additional stress it generates, whereas the stress and strain state would remain purely elastic for a matrix without inclusion. It is well known that for rolling element bearing and gear applications, the ruin of components is often linked to cracks initiated in the vicinity of large or hard inclusions that act as stress raisers. It turned out that plastic strains tend to reduce the stress generated by the contact pressure while hard heterogeneities will increase it. As plastic strain accumulation can provide the basis for fatigue damage criteria, the second half of the paper will illustrate how the method can be used to identify and rank geometrical and material parameters that influence the location and magnitude of the maximal plastic strain.