The sub-crustal stress induced by mantle convection has been traditionally computed using the Runcorn formulae of solving the Navier-Stokes problem. The main disadvantage of this method is a limited spectral resolution (up to degree 25 of spherical harmonics) due to a divergence of the spherical harmonic expression. To improve the spectral resolution, we propose a new method of computing the horizontal components of the sub-crustal stress based on utilising the stress function with a numerical differentiation. According to the proposed method, the stress function is functionally related to the gravity and crust structure models expressed in terms of spherical harmonics, instead of directly relating the stress components with partial derivatives of these spherical harmonics. The stress components are then computed from the stress function by applying a numerical differentiation. This modification increases the degree-dependent convergence domain of the asymptotically convergent series and consequently allows computing the stress components to a higher spectral resolution, which is compatible with currently available global crustal models. We further utilise the solution to the Vening Meinesz-Moritz inverse problem of isostasy in definition of the stress function. This definition facilitates a variable crustal thickness instead of assuming only a constant value adopted in the Runcorn formulae. The crustal thickness and sub-crustal stress are then determined directly from gravity data and a crustal structure model. We apply this numerical approach to compute the sub-crustal stress globally. Regional results are also presented and discussed over study areas of oceanic subduction zones, convergent continent-to-continent collision zones and hotspots. We demonstrate that the largest (in magnitude) sub-crustal stress occurs mainly along seismically active convergent tectonic plate boundaries.