The Planar Three body Problem
Example 2. Reduction of DoF using the Center of Mass

Equations of Motion

Clear["Global`*"]
x3=-(m1 x1[t]+m2 x2[t])/m3;
y3=-(m1 y1[t]+m2 y2[t])/m3;
r12=(Sqrt[(x1[t]-x2[t])^2+(y1[t]-y2[t])^2]);
r13=(Sqrt[(x1[t]-x3)^2+(y1[t]-y3)^2]);
r23=(Sqrt[(x3-x2[t])^2+(y3-y2[t])^2]);

En=(1/2)m1(px1[t]^2+py1[t]^2)+(1/2)m2(px2[t]^2+py2[t]^2)+(1/2)m3(px3^2+py3^2)-((m1 m2)/r12)-((m1 m3)/r13)-((m2 m3)/r23) ;
eq[1]=px1[t]; eq[2]=px2[t]; eq[3]=py1[t]; eq[4]=py2[t];
eq[5]=-m3((x1[t]-x3)/(r13^3))-m2((x1[t]-x2[t])/(r12^3));  (*px1*)
eq[6]=-m3((x2[t]-x3)/(r23^3))-m1((x2[t]-x1[t])/(r12^3)); (*px2*)
eq[7]=-m3((y1[t]-y3)/(r13^3))-m2((y1[t]-y2[t])/(r12^3));  (*py1*)
eq[8]=-m3((y2[t]-y3)/(r23^3))-m1((y2[t]-y1[t])/(r12^3)); (*py2*)
DifferentialEquations={x1'[t]==eq[1],x2'[t]==eq[2], y1'[t]==eq[3], y2'[t]==eq[4], 
px1'[t]==eq[5],px2'[t]==eq[6],py1'[t]==eq[7],py2'[t]==eq[8]}

Parameters and Initial Conditions

m1=1; m2=1; m3=1;
x10=0; 
x20=1;
x30=-(m1*x10+m2*x20)/m3;
y10=0.1;
y20=0;
y30=-(m1*y10+m2*y20)/m3;
px10=0;
px20=0;
px30=-(m1*px10+m2*px20)/m3;
py10=0;
py20=1;
py30=-(m1*py10+m2*py20)/m3;
tmax=20;
InitialConditions={x1[0]==x10, x2[0]==x20, y1[0]==y10, y2[0]==y20,px1[0]==px10, px2[0]==px20, py1[0]==py10, py2[0]==py20}
  Print["3rd body  x30=",x30,"  y30=",y30,"   px30=",px30,"   py30=",py30]

Numerical Solution of the system

DSystem=Union[DifferentialEquations,InitialConditions];
sol=NDSolve[DSystem,{x1,x2,y1,y2,px1,px2,py1,py2},{t,0,tmax}]
x1t[t_]=x1[t]/.sol[[1]];
x2t[t_]=x2[t]/.sol[[1]];
y1t[t_]=y1[t]/.sol[[1]];
y2t[t_]=y2[t]/.sol[[1]];
px1t[t_]=px1[t]/.sol[[1]];
px2t[t_]=px2[t]/.sol[[1]];
py1t[t_]=py1[t]/.sol[[1]];
py2t[t_]=py2[t]/.sol[[1]];
x3t[t_]=-(m1 x1t[t]+m2 x2t[t])/m3;
y3t[t_]=-(m1 y1t[t]+m2 y2t[t])/m3;
px3t[t_]=-(m1 px1t[t]+m2 px2t[t])/m3;
py3t[t_]=-(m1 py1t[t]+m2 py2t[t])/m3;

Presentation of the Numerical Solution

Energy Conservation

Ener[t_]=En/.{x1[t]->x1t[t], y1[t]->y1t[t], x2[t]->x2t[t], y2[t]->y2t[t],px1[t]->px1t[t], py1[t]->py1t[t],px2[t]->px2t[t], py2[t]->py2t[t],px3->px3t[t], py3->py3t[t]} ;
Plot[Ener[t]-Ener[0],{t,0,tmax}]

Evolution in time

Plot[{x1t[t],x2t[t],x3t[t]},{t,0,tmax},PlotStyle->{Red,Green,Blue}]
Plot[{y1t[t],y2t[t],y3t[t]},{t,0,tmax},PlotStyle->{Red,Green,Blue}]

Orbits

orbit=ParametricPlot[{{x1t[t],y1t[t]},{x2t[t],y2t[t]},{x3t[t],y3t[t]}},{t,0,tmax},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[0,0,1]}]

Clear[t0];
Animate[body1={{x1t[t],y1t[t]}} /.t->t0;
body2={{x2t[t],y2t[t]}} /.t->t0;
body3={{x3t[t],y3t[t]}} /.t->t0;
p1=ListPlot[body1, PlotRange->{{-2,2},{-2,2}}, PlotStyle->{PointSize[0.03], Red}];
p2=ListPlot[body2, PlotRange->{{-2,2},{-2,2}}, PlotStyle->{PointSize[0.03], Green}];
p3=ListPlot[body3, PlotRange->{{-2,2},{-2,2}}, PlotStyle->{PointSize[0.03], Blue}];
Show[{orbit,p1,p2,p3}], {t0,0,tmax,0.01}, AnimationRate->1]