## 1 Introduction  
The solutions for the problem of settlement of rigid rectangular footings have been proposed by a variety of authors [1; 2; 3; 4; 5; 6] with various approaches. However, the limitation of which is a lack of explanation for the failure mechanism of soil medium underneath a rigid rectangular footing. By considering the vertical displacement as a result of components due to shear strain and volumetric strain of a loaded soil element, the authors denote not only the non-linear dependence between applied load and strain components but also the approach to the bearing capacity of rigid rectangular footings. To calculate shear strain and volumetric strain, the elastic-perfectly plastic model of [PERSON] [7] and the non-linear model proposed by [PERSON] are applied to indicate the non-linear behavior of the soil subgrade below a rectangular footing subjected to a uniformly distributed load, thereby achieving a dependency of deformation on load with a double curvature. So far, the non-linear analysis of soil subgrade under a rectangular foundation has solely been conducted using the finite element method. Therefore, it is essential to consider the non-linear soil deformation properties in analytical solutions based on the results from compression and triaxial tests. As a result, the varying values of shear modulus and bulk modulus depending on the deviator stress are elucidated clearly at any point in Cartesian coordinates. In addition, linear behavior can also be observed in the non-linear solution, when considering the shear modulus and bulk modulus as constant values. Consequently, the non-linear analysis turns out to be a classic linear analysis, based on [PERSON]'s law (\(G=const\) and K = const). In this article, for the purpose of estimating the settlement and bearing capacity, we applied [PERSON]'s system of equations, which enables engineers to consider both linear and non-linear behaviors of soil media under a rectangular footing as well as to thoroughly separate the total vertical strain \(\varepsilon_{\text{z}}\) into shear strain \(\varepsilon_{\text{\tiny{Y}}}\) and volumetric strain \(\varepsilon_{\text{\tiny{V}}}\) as follows:  
\[\varepsilon_{\text{z}}=\varepsilon_{\text{\tiny{Y}}}+\varepsilon_{\text{\tiny{V} }}=\frac{\sigma_{z}-\sigma_{\text{m}}}{G\big{(}\sigma_{\text{m}},\tau_{\text{i }}/\tau_{\text{i}}^{*}\big{)}}+\frac{\sigma_{\text{m}}}{K(\sigma_{\text{m}})} \tag{1}\]  
in which \(G(\sigma_{\text{m}},\tau_{\text{i}}/\tau_{\text{i}}^{*})\) and \(K(\sigma_{\text{m}})\) are, respectively, shear modulus and bulk modulus, the values of which are relevant to the mean stress \(\sigma_{\text{m}}\) as well as the ratio of mobilized shear stress \(\tau_{\text{i}}\) to maximum shear stress \(\tau_{\text{i}}^{*}\), which can be found using the Mohr-Coulomb criterion at the considered points. In a particular case, when \(G=const\) and \(K=const\), Equation (1) turns into the Hooke's law equation [8] that is presented as:  
\[\varepsilon_{\text{z}}=\frac{\sigma_{\text{z}}-\sigma_{\text{m}}}{2G}+\frac{ \sigma_{\text{m}}}{K} \tag{2}\]  
\(G(\sigma_{\text{m}},\tau_{\text{i}}/\tau_{\text{i}}^{*})\) and \(K(\sigma_{\text{m}})\) can be retrieved from standard triaxial tests, the general view of which is observed in Figure 1 below. A significant advantage of using [PERSON]'s system of equations is that it allows for the determination of the typical characteristics of total strain, volume, and shear strain of loaded soil elements. Figure 1 accurately depicts real soil behavior, with a decaying dependency between volumetric strain and mean stress, and a persistent dependency between shear strain and deviator stress. The final state following the latter dependency leads soil elements to failure.