## 1 Introduction
One model was proposed in [10; 11; 26; 27; 28], where the isotropic loss of mass was interpreted as additional force given by  
\[\vec{F}=-\frac{\dot{m}}{2m}\dot{\vec{r}}, \tag{2}\]  
where the dot means derivation with respect to time. Considering this force as a perturbation of an aperiodic motion on a conic section, one can derive a system of differential equations determining the variation of the orbital elements. These equations were used to model the motion of planets under stellar mass loss without taking into account the mutual planetary interaction (see [29; 30; 31; 32]). From the other side, one can introduce the perturbing force in the form  
\[\vec{F}=-\frac{\dot{\gamma}(t)}{\gamma(t)}\vec{r}, \tag{3}\]  
where the function \(\gamma(t)\) is assumed to be twice differentiable and given by  
\[\gamma(t)=\frac{m_{00}+m_{10}}{m_{0}(t)+m_{1}(t)}, \tag{4}\]  
where \(m_{00}=m_{0}(t_{0})\), \(m_{10}=m_{1}(t_{0})\) are the masses of the bodies \(P_{0}\), \(P_{1}\), respectively, at the initial instant of time \(t_{0}\). In this case, the perturbed motion of the bodies may be describedin terms of the osculating elements of aperiodic motion on quasi-conic sections (see [33; 34; 35; 36]). This approach is easily generalized to the systems of three and more interacting bodies of variable masses (see [37; 38; 39; 40]). It should be emphasized that the potential energy of the many-body system is reduced to a sum of pairwise gravitational attractions of two bodies. Therefore, in a planetary system, where the mass of a parent star is usually much greater than the masses of planets, it is natural to describe the motion of the planets about the star, to the first approximation, by the corresponding solution to the two-body problem. As Equation (1) describing the relative motion of two bodies is not integrable, in general, one needs to generalize perturbation theory to cases in which the masses of the two bodies may vary independently at different rates and in which interaction with other bodies, and the reactive forces arising in the case of anisotropic mass variation, may be taken into account. The study of such models is usually associated with cumbersome symbolic calculations, for which it is expedient to use computer algebra systems, for example, Wolfram Mathematica [41]. In the present paper, we discuss a general approach to solving the problem of two bodies with variable masses based on the application of perturbation theory. We assume that the masses of the bodies change isotropically and independently of each other; therefore, the relative motion of the bodies is as described by Equation (1). Two kinds of perturbing forces, given by (3) and (4), are used to modify Equation (1) and to reduce it to a form that is appropriate for the application of perturbation theory. In the two cases, the differential equations determining the perturbed motion of the bodies are obtained in the framework of [PERSON]'s formalism in terms of the osculating elements of aperiodic motion on conic and quasi-conic sections, respectively. These equations are solved numerically, and their solutions are compared for different laws of the masses' change. For small eccentricities of the orbits, the perturbing forces are expanded in power series in small parameters up to the second-order terms. Then, the equations of the perturbed motion are averaged over the mean longitude to obtain the differential equations describing the evolution of orbital elements over long periods of time. The obtained results demonstrate that the second model, based on the aperiodic motion on a quasi-conic section, is more appropriate for generalization to the many-body problem with variable masses. To carry out all relevant symbolic and numerical calculations, we use here the computer algebra system Wolfram Mathematica [41].