Restricted TBP - evolution
v.2022

Parameters and  ODEs

Clear["Global`*"]

m=0.0123;
a=(m/3)^(1/3)   (* Hill radius *)
Body1=ListPlot[{{-m,0}},PlotStyle->{Magenta,PointSize[0.05]}];
Body2=ListPlot[{{1-m,0}},PlotStyle->{Blue,PointSize[0.02]}];

r1=(Sqrt[(x[t]+m)^2+y[t]^2]); r2=(Sqrt[(x[t]-1+m)^2+y[t]^2]);
deq1=x''[t]-2 y'[t]-x[t]==-(1-m)(x[t]+m)/r1^3-m*(x[t]-1+m)/r2^3;
deq2=y''[t]+2x'[t]-y[t]==-y[t]((1-m)/r1^3+m/r2^3);

Evolution

Initial conditions

tmax=10;
h=-1.63
x0=0.1; 
y0=0; 
vx0=0; 
vy0=3.8;
ivp={deq1,deq2,x[0]==x0, y[0]==y0,x'[0]==vx0, y'[0]==vy0}

R1=r1/.{x[t]->x0,y[t]->y0};
R2=r2/.{x[t]->x0,y[t]->y0};
Body3start=ListPlot[{{x0,y0}},PlotStyle->{Red,PointSize[0.01]}];

Numerical integration & orbit in rotating frame

sol=NDSolve[ivp,{x,y},{t,0,tmax}, MaxSteps->100000];
xx=x[t]/.sol[[1]]
yy=y[t]/.sol[[1]]
Body3end=ListPlot[{{xx/.t->tmax,yy/.t->tmax}},PlotStyle->{Gray,PointSize[0.01]}];
orbit=ParametricPlot[{xx,yy},{t,0,tmax}]

Show[{orbit, Body1,Body2, Body3start, Body3end}, PlotRange->All, AxesOrigin->{0,0}]

Orbit in Inertial Frame

XX=xx*Cos[t]-yy*Sin[t];
YY=xx*Sin[t]+yy*Cos[t];
ParametricPlot[{XX,YY},{t,0,tmax}, PlotRange->{{-2.5, 2.5},{-2.5,2.5}}]