## 3 Results  
All the equations we have derived contain the function \(\gamma(t)\) and its derivatives, which are determined by the laws of mass variation \(m_{0}(t),m_{1}(t)\). It should be emphasized that in the case of constant masses (\(\gamma(t)=1\)), all the equations coincide with the well-known results for the two-body problem of constant mass (see [3]). To find solutions to the differential Equations (17)-(22) and (29)-(34) and to compare Models 1 and 2, one needs to choose some realistic laws for the mass variation. We consider here the Eddington-Jeans law (see [24; 25]), which is used most often in applications and is given by  
\[\frac{dm}{dt}=-\alpha m^{n}, \tag{41}\]  
where \(\alpha\) is a small positive constant, and the exponent \(n\) lies between 0.4 and 4.4. Choosing \(n=2\), for example, and solving Equation (41), we obtain  
\[m(t)=\frac{m(0)}{1+\alpha m(0)t}. \tag{42}\]  
Note that parameters \(\alpha\) and \(n\) may be chosen independently for the bodies \(P_{0}\), \(P_{1}\). Furthermore, we assume that the most massive body \(P_{0}\) loses its mass according to the Eddington-Jeans law, while the mass of the body \(P_{1}\) remains constant (\(m_{1}(t)=m_{10}=const\)). As it is not possible to find analytical solutions to the differential equations under consideration, we obtain their numerical solutions using the dimensionless variables for convenience. For example, we use the initial values of the semi-major axis \(a(0)=a_{0}\) and the initial mass \(m_{00}\) of body \(P_{0}\) as units of distance and mass, respectively, and define dimensionless distance \(a_{j}^{*}\), mass \(m_{j}^{*}\), and time \(t^{*}\) by  
\[a^{*}=\frac{a}{a_{0}},\ m_{j}^{*}=\frac{m_{j}}{m_{00}},\ t^{*}=t\frac{\sqrt{ \kappa}}{a_{0}^{3/2}},\ j=0,1.\]  
The corresponding function \(\gamma(t^{*})\) takes the form  
\[\gamma(t^{*})=\frac{1+m_{10}^{*}}{1+m_{10}^{*}+\alpha^{*}t^{*}}.\]  
We also choose the following initial values for orbital elements:  
\[m_{10}^{*}=1/100,\ a^{*}(0)=1,\ e_{0}=1/20,\ i_{0}=0,\ \Omega_{0}=0,\ \omega_{0}=0,M(0)=0.\]Solutions of two systems of differential equations, Equations (17)-(22) and (29)-(34), are shown in Figures 1-4. Comparison of the semi-major axis \(a^{*}(t^{*})\) computed in the framework of Model 1 and the function \(\gamma(t^{*})\) (see Figure 1) shows that the value of \(a^{*}\) oscillates slightly near the curve \(\gamma(t^{*})\), which increases monotonically. As the secular part of the semi-major axis increases as a function of \(\gamma(t)\) (see (38)), such behavior of the function \(a^{*}(t^{*})\) seems to be correct and correspond to Equations (17) and (36). Figure 2 shows that in Model 1, the eccentricity only oscillates near its initial value \(e=1/20\), and the amplitude of these oscillations is very small; its secular part is equal to zero. Such behavior of the eccentricity corresponds to Equations (18) and (36). In Model 2, the right-hand side of the corresponding differential Equation (39) also contains oscillating terms, which are proportional to the second derivative of the function \(\gamma(t)\). As this derivative is of order \(\alpha^{2}\) (see (41)), the amplitude of oscillations is negligible, and they are not observable on Figure 2.  
Similar results are obtained for the longitude of pericenter \(\omega(t^{*})\) (see Figure 3). It should be noted only that according to the corresponding differential Equation (40), the function \(\omega(t^{*})\) should contain a non-zero secular part. However, it is of order \(\alpha^{2}\) and is negligible. The same reason results in similar behavior of the eccentric anomaly, which looks similar in both models (see Figure 4). There should be small oscillations of \(E(t^{*})\), but due to the very small amplitude, we can observe only the secular parts of the eccentric anomaly on the figure. Using the same initial conditions and the chosen function \(\gamma(t)\), we can find a numerical solution to the original Equation (1); only one Cartesian coordinate \(x(t)\) is shown in Figures 5 and 6 by the red curve.