Inclusions are the primary sites for subsurface fatigue crack initiation in bearing contacts. To understand the mechanisms of subsurface crack nucleation under contact loading, a detailed description of the stress field around these inclusions is necessary. This paper presents a new approach to computing stresses in an inhomogeneous medium where inclusions are treated as inhomogeneities in a homogeneous material matrix. The approach is based on the Discrete Element (DE) Method in which the material continuum is replaced by a set of rigid discrete interacting elements. The elements are connected to each other along their sides through springs and dampers to form the macro-continuum and undergo relative displacements in accordance with Newton’s laws of motion under the action of external loading. The spring properties are derived in terms of the overall elastic properties of the continuum. The relative motion between elements gives rise to contact forces due to stretching or compression of the inter-element springs. These forces are evaluated at each time-step and the corresponding equations of motion are solved for each element. Stresses are calculated from the inter-element joint forces. A Hertzian line contact case, with and without the presence of subsurface inclusions, is analyzed using the DE model. The DE model was used to determine stresses for an inclusion-free medium that compares well with that obtained from the continuum elasticity models. Parametric studies are then carried out to investigate the effects of size, location, orientation, and elastic properties of inclusions on the subsurface stress field. Both inclusions that are stiffer and/or softer than the base material are seen to give rise to stress concentrations. For inclusions that are stiffer than the base material (semi-infinite domain), the stress concentration effect increases with their elastic modulus. The stress concentration effect of a softer inclusion is higher than that of a stiffer inclusion. Inclusions that are oriented perpendicular to the surface give rise to much higher von Mises stresses than the ones that are oriented parallel to the surface. There is little change in the maximum von Mises stress for inclusions that are located deep within the surface.