In order to calculate the stresses, the set of equations is solved that includes the elastic equilibrium equation of the gravitating planet, and the Poisson equation that relates the anomalous gravitational potential to the density distribution. The Green’s function method (or loading numbers method) allows us to find boundary conditions based on the data on topography and the gravitational field of the planet. Relief and density anomalies distributed in special ways in the planetary interior are the sources of the anomalous gravitational field. In the stress calculation in [5], for the sake of definiteness of the problem, an assumption was made that there were two sources of the Martian anomalous gravitational field: non-equilibrium relief and density anomalies at the crust−mantle boundary (two-level compensation model). In this study, an exploration of zones of maximum shear stresses and tension−compression stresses in the Martian interior is proposed for a three-level compensation model (the non-equilibrium relief, density anomalies at the crust−mantle boundary, and density anomalies at the base of the lithosphere).