Abstract
The problem of deriving differential equations determining secular perturbations of orbital elements in a multiplanetary system is considered in the case when the central star loses its mass isotropically, while the masses of the planets can change anisotropically, which leads to the appearance of reactive forces. As a model of a multiplanetary system, the classical problem of $$n + 1$$ variable-mass bodies is used, when $$n$$ bodies move around the central star along quasi-elliptical non-intersecting orbits and interact with each other in accordance with the law of universal gravitation. It is assumed that the masses of the bodies change at different rates and the laws of mass change are considered to be arbitrary given functions of time. Differential equations of motion of bodies in osculating elements of aperiodic motion along quasi-conical orbits corresponding to the planetary Lagrange equations are derived. An algorithm for calculating the perturbing functions in the form of power series in small parameters and the derivation of differential equations determining secular perturbations of orbital elements are discussed. All necessary symbolic calculations are performed using the Wolfram Mathematica computer algebra system.