Purpose. To investigate the features of the algorithm implementation for finding the derivatives of the spatial distribution function of the planet's masses with the use of high-order Stokes constants and, on the basis of this, to find its analytical expression. According to the given methodology, to carry out calculations with the help of which to carry on the study of dynamic phenomena occurring inside an ellipsoidal planet. The proposed method involves the determination of the derivatives of the mass distribution function by the sum, the coefficients of which are obtained from the system of equations, which is incorrect. In order to solve it, an error-resistant method for calculating unknowns was used. The implementation of the construction is carried out in an iterative way, while for the initial approximation we take the three-dimensional function of the density of the Earth's masses, built according to Stokes constants up to the second order inclusive, by dynamic compression by the one-dimensional density distribution, and we determine the expansion coefficients of the derivatives of the function in the variables to the third order inclusive. They are followed by the corresponding density function, which is then taken as the initial one. The process is repeated until the specified order of approximation is reached. To obtain a stable result, we use the Cesaro summation method (method of means).. The calculations performed with the help of programs that implement the given algorithm, while the achieved high (ninth) order of obtaining the terms of the sum of calculations. The studies of the convergence of the sum of the series have been carried out, and on this basis, a conclusion has been made about the advisability of using the generalized finding of the sums based on the Cesaro method. The optimal number of contents of the sum terms has been chosen, provides convergence both for the mass distribution function and for its derivatives. Calculations of the deviations of mass distribution from the mean value ("inhomogeneities") for extreme points of the earth's geoid, which basically show the total compensation along the radius of the Earth, have been performed. For such three-dimensional distributions, calculations were performed and schematic maps were constructed according to the taken into account values of deviations of three-dimensional distributions of the mean ("inhomogeneities") at different depths reflecting the general structure of the Earth's internal structure. The presented vector diagrams of the horizontal components of the density gradient at characteristic depths (2891 km - core-mantle, 700 km - middle of the mantle, also the upper mantle - 200, 100 km) allow us to draw preliminary conclusions about the global movement of masses. At the same time, a closed loop is observed on the “core-mantle” edge, which is an analogy of a closed electric circuit. For shallower depths, differentiation of vector motions is already taking place, which gives hope for attracting these vector-grams to the study of dynamic motions inside the Earth. In fact, the vertical component (derivative with respect to the z variable) is directed towards the center of mass and confirms the main property of mass distributions - growth when approaching the center of mass. The method of stable solution of incorrect linear systems is applied, by means of which the vector-gram of the gradient of the mass distribution function is constructed. The nature of such schemes provides a tool for possible causes of mass redistribution in the middle of the planet and to identify possible factors of tectonic processes in the middle of the Earth, i.e indirectly confirms the gravitational convection of masses. The proposed technique can be used to create detailed models of density functions and its characteristics (derivatives) of the planet's interior, and the results of numerical experiments - to solve tectonics problems.