Numerics etc.

Here is Mathematica solving a real orbit problem.

Clear[f1,y,t]
f1[t_,{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=1/81.45,d1 ,d2},
d1 =((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2)); d2=((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2)) ;
{y4,y5,y6,
2y5 + y1 - ((μ(y1+μ-1))/d1) -(((1-μ)(y1+μ))/d2),
-2y4 + y2 - ((μ y2)/d1) -(((1-μ)y2)/d2),
- ((μ y3)/d1) -(((1-μ)y3)/d2)
}
]
IC = {0.994,0,0,0,-2.001585106370825224,0};
TMax=17.06521656015796; 
ySol[t_]=y[t]/.NDSolve[{
y'[t]==f1[t,y[t]],
y[0]==IC
},
y[t],
{t, 0, TMax}
][[1]];
MmaPic=ParametricPlot[ ySol[t][[{1,2}]],{t,0,TMax},
AxesLabel->{"y1","y2"}]

Euler[{t_,y_},h_,f_]:={t+h,y+h*f[t,y]}

ImprovedEuler[{t_,y_},h_,f_]:=Module[
{yStar},
yStar= y+h*f[t,y];

{t+h, y+(h/2)*(f[t,y]+f[t+h,yStar])}
]

RK4[{t_,y_},h_,f_]:=Module[
{k1,k2,k3,k4},
k1=f[t,y];
k2=f[t+(h/2),y+ (h/2)k1];
k3=f[t+(h/2),y+ (h/2)k2];
k4=f[t+h,y+ h k3];
{t+h, y+(h/6)*(k1+2k2+2k3+k4)}
]

Here are some numerical solutions

TTest=1.0;
n=5000;
h=N[TTest/n];

{t,y}={0, IC};
EulerData=Table[
{t,y}=Euler[{t,y},h,f1];
y[[{1,2}]],
{n}];

{t,y}={0,IC};
ImprovedEulerData=Table[
{t,y}=ImprovedEuler[{t,y},h,f1];
y[[{1,2}]],
{n}];

{t,y}={0,IC};
RK4Data=Table[
{t,y}=RK4[{t,y},h,f1];
y[[{1,2}]],
{n}];

Show[
MmaPic,
ListPlot[{EulerData,ImprovedEulerData,RK4Data},
PlotStyle->{Red, Blue, Green}]
]