Evaluation of stress distribution through soil strata is one of the fundamentals of geotechnical analysis. Vertical stress is attributed to the loading intensity and geostatic stress. Even though many analytical methods have been developed idealizing the soil as a homogeneous, isotropic, and elastic material, the soil consists of layers of different characteristics in reality. In this study, stress distribution under a strip footing through a two-layered system consisting of purely sandy and purely clayey soils where the sand at the top was analyzed with the finite element method using Plaxis 2D geotechnical software. The Mohr–Coulomb model was used to feed the characteristics of the soils to the computer software. Different test cases were simulated by changing the soil strength parameters; friction angle of sandy soils, undrained cohesion of clayey soils, and the thickness of two layers for selected soil types, while the total depth of the soil system was 8 m. The numerical results were compared with Boussinesq’s solutions which are linear elastic solutions. It shows that stress at a selected point is small when the soils are loose and soft compared to the soils which are dense and stiff. The lower the thickness of the sand layer higher the stress at a selected point. Vertical stress values of Boussinesq’s solutions agree with the finite element analysis values when moving away from the footing both in horizontal and vertical directions.