This research aims to develop a practical mathematical model for estimating the Equivalent Elastic Modulus (Eeq) of multi-layered soils beneath shallow foundations, overcoming the limitations of traditional experimental methods in characterizing complex soil properties. This study employs an innovative approach that combines the results of 3D Finite Element Method (FEM) modeling with the data-mining capabilities of the Evolutionary Polynomial Regression (EPR) algorithm. The research methodology is based on running 3D FEM simulations under various geometric and parametric soil conditions to generate a precise training dataset. The EPR algorithm was subsequently used to identify structured and simple relationships between the input variables and the Eeq value extracted from the FEM results. Key findings revealed that three primary factors significantly influence the equivalent modulus: the ratio of the first layer thickness to the footing width (h1/B), the ratio of the second layer thickness to the footing width (h2/B), and the elastic moduli (EEE) of the soil layers. Utilizing these non-dimensional parameters enabled the elimination of the footing width variable and the simplification of the final model. The results demonstrate that the developed EPR model can estimate Eeq with suitable accuracy for conditions corresponding to a 2.5 cm settlement at the footing center. Validation against independent data proved the high generalizability of the derived relationship. The final conclusion is that the EPR model, by transforming complex numerical computations into a simple and practical mathematical relationship, can serve as an effective tool for geotechnical engineers to reduce routine calculations associated with shallow foundations, especially in scenarios requiring accurate estimation of multi-layered soil properties.