Sub-lithospheric stress due to mantle convection can be determined from gravimetric data based on Runcorn’s theory. In this paper, the satellite gradiometric data of the recent European satellite mission, the Gravity field and steady-state Ocean Circulation Explorer (GOCE) is used to determine the sub- lithospheric stress locally in Iran. The method of S function (SF) with numerical differentiation is developed further and an integral equation connecting satellite gradiometric data to SF is presented. The integral equation will be used to invert the real gradiometric data of GOCE to recover the SF. Later on, the sub-lithospheric shear stresses, which are the northwards and eastwards derivatives of the SF, are computed numerically. Our numerical results show that the mean square error of the recovered SF is smaller than the values of the SF meaning that the recovery process is successful. Also, the recovered stress has a good agreement with the tectonic boundaries and active seismic points of the world stress map (WSM) database. This stress reaches amplitude of 100 MPa in the territory.