Проверь, работает ли этот код
c0 = 299792458.; G = 6.67384*10^(-11); GM =
G*1.9891*10^(30); a = 1950100000; T1 = 7.7519387738647*3600; t0 =
T1*Sqrt[GM/a^3]; e0 = 0.6171334; m2 = 1.387; m1 = 1.441; v0 =
450000/Sqrt[GM/a]*2.003857; c =
c0/Sqrt[GM/a]; x10 = (1 - e0)*m2/(m1 + m2); x20 = -(1 - e0)*
m1/(m1 + m2); y10 = 0.; y20 = 0; v1y =
v0*m2/(m1 + m2); v2y = -v0*m1/(m1 + m2); v1x = 0.; v2x = 0;
r1 = {x1[t], y1[t], z1[t]};
r2 = {x2[t], y2[t], z2[t]};
v1 = D[r1, t];
v2 = D[r2, t];
r = r1 - r2;
ra = Norm[r];
rn = r/ra;
eq1 = -D[v1, t] -
m2*rn/ra^2 + (m2/(c*ra)^2)*(rn*((5*m1 + 4*m2)/ra - v1.v1 -
2*v2.v2 + 4*v1.v2 + (3/2)*(v2.rn)^2) + (v1 -
v2)*(rn.(4*v1 - 3*v2)));
eq2 = -D[v2, t] +
m1*rn/ra^2 - (m1/(c*ra)^2)*(rn*((5*m2 + 4*m1)/ra - v2.v2 -
2*v1.v1 + 4*v1.v2 + (3/2)*(v1.rn)^2) - (v1 -
v2)*(rn.(4*v2 - 3*v1)));
eq0 = {x1[0] == x10, y1[0] == 0, z1[0] == 0, x1'[0] == 0,
y1'[0] == v1y, z1'[0] == 0, x2[0] == x20, y2[0] == 0, z2[0] == 0,
x2'[0] == 0, y2'[0] == v2y, z2'[0] == 0};
s = NDSolve[{Table[eq1[[k]] == 0, {k, 1, 3}],
Table[eq2[[k]] == 0, {k, 1, 3}], eq0}, {x1[t], y1[t], z1[t],
x2[t], y2[t], z2[t], x1'[t], y1'[t], x2'[t], y2'[t]}, {t, 0, t0}];
{ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}} /. First[s], {t, .0,
t0}, AxesLabel -> {"x", "y"}],
ParametricPlot[{x1[t] - x2[t], y1[t] - y2[t]} /. First[s], {t, .0,
t0}, AxesLabel -> {"x", "y"}]}