ep = {0, 0};
me = 5;
mp = {{0, 5}, {0, 5}, {0, 5}};
acc = {me (ep - mp[[1]])/(Norm[ep - mp[[1]]]^1 Norm[ep - mp[[1]]]^2), me (ep - mp[[2]])/(Norm[ep - mp[[2]]]^1 Norm[ep - mp[[2]]]^1.9), me (ep - mp[[3]])/(Norm[ep - mp[[3]]]^1 Norm[ep - mp[[3]]]^2.1)};
mv = mv = {{Sqrt[Abs[Norm[acc[[1]] ] Norm[mp[[1]] - ep] ]], 0.2}, {Sqrt[Abs[Norm[acc[[1]] ] Norm[mp[[1]] - ep] ]], 0.2}, {Sqrt[Abs[Norm[acc[[1]] ] Norm[mp[[1]] - ep] ]], 0.2}};
dt = 0.2;
mpold = mp;
mp = mpold + mv dt + acc/2 dt^2;
evo = Reap[Do[
acc = {me (ep - mp[[1]])/(Norm[ep - mp[[1]]]^1 Norm[ep - mp[[1]]]^2),
me (ep - mp[[2]])/(Norm[ep - mp[[2]]]^1 Norm[ep - mp[[2]]]^1.9),
me (ep - mp[[3]])/(Norm[ep - mp[[3]]]^1 Norm[ep - mp[[3]]]^2.1)};
mpoldold = mpold;
mpold = mp;
mp = 2 mpold - mpoldold + acc dt^2;
Sow[mp];
, {1500}];][[2, 1]];
plots = Table[
Legended[
Graphics[{Gray, Disk[ep, 0.1 ],
Purple, Disk[evo[[j, 2]], 0.5 ], Line[evo[[1 ;; j, 2]] ]
,
Orange, Disk[evo[[j, 3]], 0.5 ], Line[evo[[1 ;; j, 3]] ]
,
Black, Disk[evo[[j, 1]], 0.5 ], Line[evo[[1 ;; j, 1]] ]
},
PlotRange -> {{-10, 10}, {-10, 10}}, Frame -> False], LineLegend[{Black, Purple, Orange}, {"F\[Proportional]\!\(\*FractionBox[\(1\), \SuperscriptBox[\(r\), \(2\)]]\)", "F\[Proportional]\!\(\*FractionBox[\(1\), SuperscriptBox[\(r\), \(1.9\)]]\)", "F\[Proportional]\!\(\*FractionBox[\(1\), SuperscriptBox[\(r\), \(2.1\)]]\)"}] ]
, {j, 1, Dimensions[evo][[1]]}];
ListAnimate[plots[[1 ;; -1 ;; 5]] ]