Linearization of the nonlinear equations and iterative solution is the most well-known scheme in many engineering problems. For geodetic applications of the LEO satellites, e.g. the Earth’s gravity field recovery, one needs to provide an initial guess of the satellite location or the so-called reference orbit. Numerical integration can be utilized for generating the reference orbit if a satellite’s state vector i.e., position and velocity is known at a reference epoch. However, the numerically integrated orbit deviated from the real orbit due to imperfect force models. The more accurate reference orbit the less linearization error occurs. The deviation between the reference and real orbit can be minimized using the least squares method. Different analytical and numerical techniques have been developed for calculation of the design matrix of the least squares method. Herein, we have generalized the idea of the Lagrange coefficients for the determination of the design matrix’s entries in the gravitational field an attracting inhomogeneous mass body. Numerical implementation of the proposed method shows its high performance.