A numerical study of the mechanical behavior of heterogeneous soil systems, consisting of a bulk of sand with embedded stiff gravel inclusions or soft clay inclusions, is performed. A solution scheme using parallel computing is employed when analyzing two different categories of problems. First, a homogenization problem is studied, where use of a single representative volume element subjected to plane strain compression offers the possibility to investigate the coupling between the response at a local scale and at a global scale. Second, a plane strain footing problem with different heterogeneous soil systems is analyzed using a traditional finite element formulation. The material model utilized for the soil is a large deformation formulation of non-associated elasto-plasticity with an isotropic hardening law, able to represent dilation. It was found that the shape of the gravel or clay inclusions in the systems had no significant effect on the global responses, whereas the strain localizations in the two different soil systems, sand-gravel and sand-clay, were found to have different character. The effect of the initial density on the response was clearly observed in the localization patterns.