(*1. 各惑星の微分方程式を数値的に解く*)sol=NDSolve[{(*惑星1（水色：地球のような円軌道に近い惑星）*)x1''[t]==-x1[t]/(x1[t]^2+y1[t]^2)^(3/2),y1''[t]==-y1[t]/(x1[t]^2+y1[t]^2)^(3/2),x1[0]==1,y1[0]==0,x1'[0]==0,y1'[0]==1,(*惑星2（オレンジ：火星のような少し外側を回る惑星）*)x2''[t]==-x2[t]/(x2[t]^2+y2[t]^2)^(3/2),y2''[t]==-y2[t]/(x2[t]^2+y2[t]^2)^(3/2),x2[0]==1.5,y2[0]==0,x2'[0]==0,y2'[0]==0.8,(*天体3（白：ハレー彗星のような細長い楕円軌道）*)x3''[t]==-x3[t]/(x3[t]^2+y3[t]^2)^(3/2),y3''[t]==-y3[t]/(x3[t]^2+y3[t]^2)^(3/2),x3[0]==2,y3[0]==0,x3'[0]==0,y3'[0]==0.4},{x1,y1,x2,y2,x3,y3},{t,0,50}];

(*2. アニメーションの作成*)
Animate[ParametricPlot[{{Evaluate[x1[s]/. sol],Evaluate[y1[s]/. sol]},{Evaluate[x2[s]/. sol],Evaluate[y2[s]/. sol]},{Evaluate[x3[s]/. sol],Evaluate[y3[s]/. sol]}},{s,0,t+0.0001},PlotRange->{{-2.5,2.5},{-2.5,2.5}},PlotStyle->{{LightBlue,Dashed},{Orange,Dashed},{White,Dotted}},Background->Black,Axes->False,(*宇宙空間っぽくするために軸を消去*)Epilog->{(*太陽*)Yellow,PointSize[0.08],Point[{0,0}],(*惑星1*)LightBlue,PointSize[0.04],Point[{Evaluate[x1[t]/. sol][[1]],Evaluate[y1[t]/. sol][[1]]}],(*惑星2*)Orange,PointSize[0.03],Point[{Evaluate[x2[t]/. sol][[1]],Evaluate[y2[t]/. sol][[1]]}],(*彗星*)White,PointSize[0.02],Point[{Evaluate[x3[t]/. sol][[1]],Evaluate[y3[t]/. sol][[1]]}]}],{t,0,48,0.1},AnimationRunning->False]