(Hier sind die Befehle für die Animationen)

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

ruku[f_,y0_,x_,nn_]:=N[NestList[rukustep[f,#,x/nn]&,y0,nn]];
ruku2[f_,y0_,x_,nn_]:=N[NestList[rukustep[f,#,(1/1000)(((#[[1,1]].#[[1,1]] +#[[1,2]].#[[1,2]]+ #[[1,3]].#[[1,3]]+#[[1,4]].#[[1,4]] +#[[1,5]].#[[1,5]]+ #[[1,6]].#[[1,6]]))^x)]&,y0,nn]];
f[{θ_,x_}]:={{θ[[4]],θ[[5]],θ[[6]],-G((m2(θ[[1]]-θ[[2]]))/((((θ[[1]]-θ[[2]]).(θ[[1]]-θ[[2]])))^(3/2)))-G((m3(θ[[1]]-θ[[3]]))/((((θ[[1]]-θ[[3]]).(θ[[1]]-θ[[3]])))^(3/2))),-G((m3(θ[[2]]-θ[[3]]))/((((θ[[2]]-θ[[3]]).(θ[[2]]-θ[[3]])))^(3/2)))-G((m1(θ[[2]]-θ[[1]]))/((((θ[[2]]-θ[[1]]).(θ[[2]]-θ[[1]])))^(3/2))),-G((m1(θ[[3]]-θ[[1]]))/((((θ[[3]]-θ[[1]]).(θ[[3]]-θ[[1]])))^(3/2)))-G((m2(θ[[3]]-θ[[2]]))/((((θ[[3]]-θ[[2]]).(θ[[3]]-θ[[2]])))^(3/2)))},1.};
f2[{θ_,x_}]:={{θ[[5]],θ[[6]],θ[[7]],θ[[8]],-G((m2(θ[[1]]-θ[[2]]))/((((θ[[1]]-θ[[2]]).(θ[[1]]-θ[[2]])))^(3/2)))-G((m3(θ[[1]]-θ[[3]]))/((((θ[[1]]-θ[[3]]).(θ[[1]]-θ[[3]])))^(3/2)))-G((m4(θ[[1]]-θ[[4]]))/((((θ[[1]]-θ[[4]]).(θ[[1]]-θ[[4]])))^(3/2))),-G((m3(θ[[2]]-θ[[3]]))/((((θ[[2]]-θ[[3]]).(θ[[2]]-θ[[3]])))^(3/2)))-G((m1(θ[[2]]-θ[[1]]))/((((θ[[2]]-θ[[1]]).(θ[[2]]-θ[[1]])))^(3/2)))-G((m4(θ[[2]]-θ[[4]]))/((((θ[[2]]-θ[[4]]).(θ[[2]]-θ[[4]])))^(3/2))),-G((m4(θ[[3]]-θ[[4]]))/((((θ[[3]]-θ[[4]]).(θ[[3]]-θ[[4]])))^(3/2)))-G((m2(θ[[3]]-θ[[2]]))/((((θ[[3]]-θ[[2]]).(θ[[3]]-θ[[2]])))^(3/2)))-G((m1(θ[[3]]-θ[[1]]))/((((θ[[3]]-θ[[1]]).(θ[[3]]-θ[[1]])))^(3/2))),-G((m3(θ[[4]]-θ[[3]]))/((((θ[[4]]-θ[[3]]).(θ[[4]]-θ[[3]])))^(3/2)))-G((m2(θ[[4]]-θ[[2]]))/((((θ[[4]]-θ[[2]]).(θ[[4]]-θ[[2]])))^(3/2)))-G((m1(θ[[4]]-θ[[1]]))/((((θ[[4]]-θ[[1]]).(θ[[4]]-θ[[1]])))^(3/2)))},1.};
f3[{θ_,x_}]:={{θ[[6]],θ[[7]],θ[[8]],θ[[9]],θ[[10]],-G((m2(θ[[1]]-θ[[2]]))/((((θ[[1]]-θ[[2]]).(θ[[1]]-θ[[2]])))^(3/2)))-G((m3(θ[[1]]-θ[[3]]))/((((θ[[1]]-θ[[3]]).(θ[[1]]-θ[[3]])))^(3/2)))-G((m4(θ[[1]]-θ[[4]]))/((((θ[[1]]-θ[[4]]).(θ[[1]]-θ[[4]])))^(3/2)))-G((m5(θ[[1]]-θ[[5]]))/((((θ[[1]]-θ[[5]]).(θ[[1]]-θ[[5]])))^(3/2))),-G((m3(θ[[2]]-θ[[3]]))/((((θ[[2]]-θ[[3]]).(θ[[2]]-θ[[3]])))^(3/2)))-G((m1(θ[[2]]-θ[[1]]))/((((θ[[2]]-θ[[1]]).(θ[[2]]-θ[[1]])))^(3/2)))-G((m4(θ[[2]]-θ[[4]]))/((((θ[[2]]-θ[[4]]).(θ[[2]]-θ[[4]])))^(3/2)))-G((m5(θ[[2]]-θ[[5]]))/((((θ[[2]]-θ[[5]]).(θ[[2]]-θ[[5]])))^(3/2))),-G((m4(θ[[3]]-θ[[4]]))/((((θ[[3]]-θ[[4]]).(θ[[3]]-θ[[4]])))^(3/2)))-G((m2(θ[[3]]-θ[[2]]))/((((θ[[3]]-θ[[2]]).(θ[[3]]-θ[[2]])))^(3/2)))-G((m1(θ[[3]]-θ[[1]]))/((((θ[[3]]-θ[[1]]).(θ[[3]]-θ[[1]])))^(3/2)))-G((m5(θ[[3]]-θ[[5]]))/((((θ[[3]]-θ[[5]]).(θ[[3]]-θ[[5]])))^(3/2))),-G((m3(θ[[4]]-θ[[3]]))/((((θ[[4]]-θ[[3]]).(θ[[4]]-θ[[3]])))^(3/2)))-G((m2(θ[[4]]-θ[[2]]))/((((θ[[4]]-θ[[2]]).(θ[[4]]-θ[[2]])))^(3/2)))-G((m1(θ[[4]]-θ[[1]]))/((((θ[[4]]-θ[[1]]).(θ[[4]]-θ[[1]])))^(3/2)))-G((m5(θ[[4]]-θ[[5]]))/((((θ[[4]]-θ[[5]]).(θ[[4]]-θ[[5]])))^(3/2))),-G((m3(θ[[5]]-θ[[3]]))/((((θ[[5]]-θ[[3]]).(θ[[5]]-θ[[3]])))^(3/2)))-G((m2(θ[[5]]-θ[[2]]))/((((θ[[5]]-θ[[2]]).(θ[[5]]-θ[[2]])))^(3/2)))-G((m1(θ[[5]]-θ[[1]]))/((((θ[[5]]-θ[[1]]).(θ[[5]]-θ[[1]])))^(3/2)))-G((m4(θ[[5]]-θ[[4]]))/((((θ[[4]]-θ[[5]]).(θ[[4]]-θ[[5]])))^(3/2)))},1.};
(*Die zu lösende Differentialgleichung*)
sol[θ0_,nn_,End_]:={#[[2]],#[[1]]}&/@ruku[f,{θ0,0.},End,nn];
sol2[θ0_,nn_,End_]:={#[[2]],#[[1]]}&/@ruku[f2,{θ0,0.},End,nn];
sol3[θ0_,nn_,End_]:={#[[2]],#[[1]]}&/@ruku[f3,{θ0,0.},End,nn];
(*Lösung der Differentialgleichung*)
KinetischeEnergie[θ_]:=(1/2)m1 θ[[4]].θ[[4]]+(1/2)m2 θ[[5]].θ[[5]]+(1/2)m3 θ[[6]].θ[[6]]
SchwerPunkt[θ_]:=(1/(m1+m2+m3))(m1 θ[[1]]+m2 θ[[2]]+m3 θ[[3]])
cr[{x1_,y1_},{x2_,y2_}]:=x1 y2-x2 y1
DrehImpuls[θ_]:=m1 cr[θ[[1]],θ[[4]]]+m2 cr[θ[[2]],θ[[5]]]+m3 cr[θ[[3]],θ[[6]]]
GesamtImpuls[θ_]:=m1 θ[[4]]+m2 θ[[5]]+m3 θ[[6]]
KinetischeEnergie2[θ_]:=(1/2)m1 θ[[5]].θ[[5]]+(1/2)m2 θ[[6]].θ[[6]]+(1/2)m3 θ[[7]].θ[[7]]+(1/2)m4 θ[[8]].θ[[8]]
SchwerPunkt2[θ_]:=(1/(m1+m2+m3+m4))(m1 θ[[1]]+m2 θ[[2]]+m3 θ[[3]]+m4 θ[[4]])
cr[{x1_,y1_},{x2_,y2_}]:=x1 y2-x2 y1
DrehImpuls2[θ_]:=m1 cr[θ[[1]],θ[[5]]]+m2 cr[θ[[2]],θ[[6]]]+m3 cr[θ[[3]],θ[[7]]]+m4 cr[θ[[4]],θ[[8]]];
GesamtImpuls2[θ_]:=m1 θ[[5]]+m2 θ[[6]]+m3 θ[[7]]+m4 θ[[8]];
Dreh[α_]:={{Cos[α],-Sin[α]},{Sin[α],Cos[α]}}

GesamtImpuls3[θ_]:=m1 θ[[6]]+m2 θ[[7]]+m3 θ[[8]]+m4 θ[[9]]+m5 θ[[10]]
SchwerPunkt3[θ_]:=(1/(m1+m2+m3+m4+m5))(m1 θ[[1]]+m2 θ[[2]]+m3 θ[[3]]+m4 θ[[4]]+m5 θ[[5]])

m1=1;
m2=1;
m3=1;
n=3000;
v=(Sqrt[((Sqrt[3])/3.)]);
s7=sol[{{1,0},Dreh[(2/3)π].{1,0},Dreh[(4/3)π].{1,0},{0,v},Dreh[(2/3)π].{0,v},Dreh[(4/3)π].{0,v}},n,30];
s9=sol[{{1+(10^(-5)),0},Dreh[(2/3)π].{1,0},Dreh[(4/3)π].{1,0},{0,v},Dreh[(2/3)π].{0,v},Dreh[(4/3)π].{0,v}},n,30];
list7={Table[s7[[i,2,1]],{i,n+1}],Table[s7[[i,2,2]],{i,n+1}],Table[s7[[i,2,3]],{i,1,n+1}]};
list9={Table[s9[[i,2,1]],{i,n+1}],Table[s9[[i,2,2]],{i,n+1}],Table[s9[[i,2,3]],{i,1,n+1}]};

listplot2=Table[0,{i,1,10}];

listplot2[[1]]=Manipulate[Show[ListPlot[list7,PlotRange->All,AspectRatio ->Automatic],Graphics[{PointSize[Large],Point[{s7[[i,2,1]],s7[[i,2,2]],s7[[i,2,3]]}]}],Axes->False],{i,1,n+1,1}];

m1=1;
m2=1;
m3=1;
n=3000;
v=0.5;
s8=sol[{{0.97000436,-0.24308753},-{0.97000436,-0.24308753},{0.,0.},-(1/2){-0.93240737,-0.86473146},-(1/2){-0.93240737,-0.86473146},{-0.93240737,-0.86473146}},n,20];
s10=sol[{{0.97000436+5 (10^(-2)),-0.24308753},-{0.97000436,-0.24308753},{0.,0.},-(1/2){-0.93240737,-0.86473146},-(1/2){-0.93240737,-0.86473146},{-0.93240737,-0.86473146}},n,20];
list8={Table[s8[[i,2,1]],{i,n+1}],Table[s8[[i,2,2]],{i,n+1}],Table[s8[[i,2,3]],{i,1,n+1}]};
list10={Table[s10[[i,2,1]],{i,n+1}],Table[s10[[i,2,2]],{i,n+1}],Table[s10[[i,2,3]],{i,1,n+1}]};

listplot2[[2]]=Manipulate[Show[ListPlot[list8,PlotRange->All,AspectRatio ->Automatic],Graphics[{PointSize[Large],Point[{s8[[i,2,1]],s8[[i,2,2]],s8[[i,2,3]]}]}],Axes->False],{i,1,n+1,1}];

m1=1;
m2=1;
m3=1;
m4=1;
m5=1;
n=3000;
v=0.5;
s11=sol3[{{1.657666,0},{0.439775,-0.169717},{-1.268608,-0.267651},{-1.268608,0.267651},{0.439775,0.169717},{0,-0.593786},{1.822785,0.128248},{1.271564,0.168645},{-1.271564,0.168645},{-1.822785,0.128248}},n,2];
list11={Table[s11[[i,2,1]],{i,n+1}],Table[s11[[i,2,2]],{i,n+1}],Table[s11[[i,2,3]],{i,1,n+1}],Table[s11[[i,2,4]],{i,1,n+1}],Table[s11[[i,2,5]],{i,1,n+1,5}]};
s2=sol2[{{ 1.382857,0},{- 1.382857,0},{0.,0.157030 },-{0.,0.157030 },{0,0.584873 },-{0,0.584873 },{1.871935,0},-{1.871935,0}},n,2];
list2={Table[s2[[i,2,1]],{i,n+1}],Table[s2[[i,2,2]],{i,n+1}],Table[s2[[i,2,3]],{i,1,n+1}],Table[s2[[i,2,4]],{i,1,n+1,5}]};

listplot2[[3]]=Manipulate[
Grid[
{
{
Show[ListPlot[list2,Axes->{False,False},PlotRange->All,AspectRatio ->Automatic],
Graphics[{PointSize[Large],Point[{s2[[i,2,1]],s2[[i,2,2]],s2[[i,2,3]],s2[[i,2,4]]}]}]]},{
Show[ListPlot[list11,Axes->{False,False},PlotRange->All,AspectRatio ->Automatic],
Graphics[{PointSize[Large],Point[{s11[[i,2,1]],s11[[i,2,2]],s11[[i,2,3]],s11[[i,2,4]],s11[[i,2,5]]}]}]]
}
}
],{i,1,n+1,1}
]

m1=1;
m2=1;
m3=(10^(-5));
m4=(10^(-6));
n=5000;
v=0.5;
s1=sol2[{{1,0},{-1,0},{0.8,0.},{7,0.},{0.,v},{0.,-v},{0.,3},{0.,0.535}},n,80.5];
list1={Table[s1[[i,2,1]],{i,n+1}],Table[s1[[i,2,2]],{i,n+1}],Table[s1[[i,2,3]],{i,1,n+1}],Table[s1[[i,2,4]],{i,1,n+1}]};

listplot2[[4]]=Manipulate[Show[Plot[{},{x,-10,10},PlotRange->{{-8,8},{-8,8}},ClippingStyle->None,AspectRatio -> 1,Axes->{False,False}],Graphics[{PointSize[Large],Yellow,Point[{s1[[i,2,1]],s1[[i,2,2]]}]}],Graphics[{PointSize[Medium],Point[{s1[[i,2,3]],s1[[i,2,4]]}]}]],{i,1,n+1,1}]

Das Notebook muss, damit die Animationen korrekt dargesetellt werden, ausgeführt werden: (am einfachsten mit EvaluateNotebook)

Das 3-Körper-Problem:

"Kreisbahn"

Die 3 Planeten haben gleiche Massen

listplot2[[1]]

"Die Acht"

Die 3 Planeten haben gleiche Massen

listplot2[[2]]

p- und s-Laufbahn

Die beiden “Sonnen” haben sehr große Masse, 
die Masse der Kometen ist vernachlässigbar klein

listplot2[[4]]

Das 4- und 5-Körper-Problem:

"Die Superacht"

Alle “Planten” haben die gleiche Masse

listplot2[[3]]

(Hier sind die Befehle für die Animationen)

listplot=Table[0,{i,9}];
Startwerte=Table[0,{i,9}];
listplot[[1]]=;
Startwerte[[1]]={{0.19991472060830104`},{0.`},{-0.005839904460257763`}};
listplot[[2]]=;
Startwerte[[2]]={{0.19991472060830104`},{0.`},{-0.005839904460257763`}};
listplot[[3]]=;
Startwerte[[3]]={{0.168294},{0.},{0.10806}};
listplot[[4]]=;
Startwerte[[4]]={{-0.18004},{0.025664},{-0.0832294}};
listplot[[5]]=;
Startwerte[[5]]={{-0.197914},{0.028212},{-0.0058399}};
listplot[[6]]=;
Startwerte[[6]]={{-0.451843},{0.213599},{-0.0145998}};
listplot[[7]]=;
Startwerte[[7]]={{-0.19282},{0.0274859},{-0.0454404}};
listplot[[8]]=;
Startwerte[[8]]={{0.186408},{0.},{0.0724716}};listplot[[9]]=;
Startwerte[[9]]={{0.0949616},{0.0401491},{-0.171378}};

a=1;
ω=1;
G =1;
M=1;
m=1;
μ=(1/(m+M)); (*Defniere die ganzen Konstanten*)
Sonne1={-μ a,0,0}; (*Ort der Sonne, für den Plot*)
Sonne2={(1-μ)a,0,0}; (*Ort der Erde, für den Plot*)
V[ξ1_,ξ2_,ξ3_]:=((-G M m(1-μ))/(Sqrt[(((ξ1+a μ))^2)+((ξ2)^2)+((ξ3)^2)]))+((-G M m μ)/(Sqrt[(((ξ1-a(1-μ)))^2)+((ξ2)^2)+((ξ3)^2)]))(*gegebenes Potential*)

Laufbahnmrs[Lb_,T_]:={
 sol=NDSolve[{(({{ξ1''[t]}, {ξ2''[t]}, {ξ3''[t]}}))==(({{2m ω ξ2'[t]+m (ω^2)ξ1[t]-D[V[ξ1[t],ξ2[t],ξ3[t]],ξ1[t]]}, {-2m ω ξ1'[t]+m (ω^2)ξ2[t]-D[V[ξ1[t],ξ2[t],ξ3[t]],ξ2[t]]}, {-D[V[ξ1[t],ξ2[t],ξ3[t]],ξ3[t]]}})),
                               (({{ξ1[0]}, {ξ2[0]}, {ξ3[0]}}))==(({{0}, {0}, {0}})),(({{ξ1'[0]}, {ξ2'[0]}, {ξ3'[0]}}))==Startwerte[[Lb]]} ,{ξ1,ξ2,ξ3},{t,0,T}];
ψ=Table[(({{ξ1[t]}, {ξ2[t]}, {ξ3[t]}})),{t,0,T,0.005}]/.sol;
Manipulate[Show[
ParametricPlot3D[
{Sonne1.(({{Cos[0.005*time], -Sin[0.005*time], 0}, {Sin[0.005*time], Cos[0.005*time], 0}, {0, 0, 1}})),Sonne2.(({{Cos[0.005*time], -Sin[0.005*time], 0}, {Sin[0.005*time], Cos[0.005*time], 0}, {0, 0, 1}})),Flatten[ψ[[1,time]]].(({{Cos[0.005*time], -Sin[0.005*time], 0}, {Sin[0.005*time], Cos[0.005*time], 0}, {0, 0, 1}}))},{t,0,T},PlotStyle->{{Thickness[0.04],Yellow},{Thickness[0.04],Yellow},{Thickness[0.01],Blue}},PlotRange->{{-1,1},{-1,1},{-0.2,0.2}}]],{time,1,(T/0.005),1}]}[[1]]

Das 3-Körper-Problem:

"exotische Trajektorien"

9 Verschiedene Laufbahnen: im Mitrotierendem-Schwerpunktskoordinaten-System
Die beiden “Sonnen” haben die gleichen Massen,
der Komet hat eine  vernachlässigbar kleine Masse

Manipulate[MatrixForm[Startwerte[[i]]];Show[listplot[[i]]],{i,1,9,1}]

Betrachte nun die erste Laufbahn, in nicht mit rotierenden Koordinaten

Laufbahnmrs[1,27.83]