(* En este primer ejemplo pueden ver tres cuerpos de la misma masa en una situación de muy alta simetría. Calculé que la fuerza centrífuga fuera compenzada exactamente por la fuerza gravitacional. Por este motivo las particulas se mueven en círculo. Los tres estan en el plano z=0. Cada una de las particulas se encuentra en un principio en uno de los vértices de un triangulo equilatero circunscrito en un circulo y todas tienen velocidad dada por v0=(Sqrt[((G m1 2 Cos[π/6])/(3 r0))]).
 Esta velocidad garantiza que la fuerza gravitacional y la centrífuga se cancelen *)

F[r_,m_,i_,j_,k_]:=-G m[[j]]m[[k]](((r[[j,i]]-r[[k,i]]))/((((r[[j]]-r[[k]]).(r[[j]]-r[[k]])))^(3/2)));
n=3;
m={m1,m2,m3};
r1={x1[t],y1[t],z1[t]};
r2={x2[t],y2[t],z2[t]};
r3={x3[t],y3[t],z3[t]};
r={r1,r2,r3};
v0=(Sqrt[((G m1 2 Cos[π/6])/(3 r0))]);
vars=Flatten[Table[r[[i]],{i,1,n}]]
syseq=Flatten[Table[m[[i]]D[r[[i,j]],{t,2}]==Sum[If[k≠i,F[r,m,j,i,k],0],{k,1,n}],{i,1,n},{j,1,3}]];
syseq=Join[syseq,{x1[0]== 0.0,y1[0]==1.0,z1[0]==0.0,x2[0]==Sin[π/3],y2[0]==-Cos[π/3],z2[0]==0,x3[0]==-Sin[π/3],y3[0]==-Cos[π/3],z3[0]==0.0,x1'[0]==-v0,y1'[0]==0.0,z1'[0]==0.0,x2'[0]==Cos[π/3]v0,y2'[0]==Sin[π/3]v0,z2'[0]==0,x3'[0]==Cos[π/3]v0,y3'[0]==-Sin[π/3]v0,z3'[0]==0}];

numsyseq=syseq/.{G->1,m1->1,m2->1,m3->1,r0->1}
sol=NDSolve[numsyseq,vars,{t,0,2.0}]
ParametricPlot3D[{(r[[1]]/.sol[[1]]),(r[[2]]/.sol[[1]]),(r[[3]]/.sol[[1]])},{t,0.0,2.0},PlotRange->All]

(* En este caso tomo el ejemplo anterior y saco a dos de las partículas fuera del
plano z=0. También dupliqué una de las masas. Como se rompieron varias de las simetrías del ejemplo anterior empiezan a pasar cosas más interesantes. *)
F[r_,m_,i_,j_,k_]:=-G m[[j]]m[[k]](((r[[j,i]]-r[[k,i]]))/((((r[[j]]-r[[k]]).(r[[j]]-r[[k]])))^(3/2)));
n=3;
m={m1,m2,m3};
r1={x1[t],y1[t],z1[t]};
r2={x2[t],y2[t],z2[t]};
r3={x3[t],y3[t],z3[t]};
r={r1,r2,r3};
v0=(Sqrt[((G m1 2 Cos[π/6])/(3 r0))]);
vars=Flatten[Table[r[[i]],{i,1,n}]]
syseq=Flatten[Table[m[[i]]D[r[[i,j]],{t,2}]==Sum[If[k≠i,F[r,m,j,i,k],0],{k,1,n}],{i,1,n},{j,1,3}]];
syseq=Join[syseq,{x1[0]== 0.0,y1[0]==1.0,z1[0]==0.0,x2[0]==Sin[π/3],y2[0]==-Cos[π/3],z2[0]==0.1,x3[0]==-Sin[π/3],y3[0]==-Cos[π/3],z3[0]==-0.2,x1'[0]==-v0,y1'[0]==0.0,z1'[0]==0.0,x2'[0]==Cos[π/3]v0,y2'[0]==Sin[π/3]v0,z2'[0]==0,x3'[0]==Cos[π/3]v0,y3'[0]==-Sin[π/3]v0,z3'[0]==0}];
numsyseq=syseq/.{G->1,m1->2,m2->1,m3->1,r0->1}
sol=NDSolve[numsyseq,vars,{t,0,6.0}]
ParametricPlot3D[{(r[[1]]/.sol[[1]]),(r[[2]]/.sol[[1]]),(r[[3]]/.sol[[1]])},{t,0.0,6.0},PlotRange->All]

(* Este es el caso anterior pero durante más tiempo. ¿Porqué crees que el sistema se mueve hacia la derecha? *)
F[r_,m_,i_,j_,k_]:=-G m[[j]]m[[k]](((r[[j,i]]-r[[k,i]]))/((((r[[j]]-r[[k]]).(r[[j]]-r[[k]])))^(3/2)));
n=3;
m={m1,m2,m3};
r1={x1[t],y1[t],z1[t]};
r2={x2[t],y2[t],z2[t]};
r3={x3[t],y3[t],z3[t]};
r={r1,r2,r3};
v0=(Sqrt[((G m1 2 Cos[π/6])/(3 r0))]);
vars=Flatten[Table[r[[i]],{i,1,n}]]
syseq=Flatten[Table[m[[i]]D[r[[i,j]],{t,2}]==Sum[If[k≠i,F[r,m,j,i,k],0],{k,1,n}],{i,1,n},{j,1,3}]];
syseq=Join[syseq,{x1[0]== 0.0,y1[0]==1.0,z1[0]==0.0,x2[0]==Sin[π/3],y2[0]==-Cos[π/3],z2[0]==0.1,x3[0]==-Sin[π/3],y3[0]==-Cos[π/3],z3[0]==-0.2,x1'[0]==-v0,y1'[0]==0.0,z1'[0]==0.0,x2'[0]==Cos[π/3]v0,y2'[0]==Sin[π/3]v0,z2'[0]==0,x3'[0]==Cos[π/3]v0,y3'[0]==-Sin[π/3]v0,z3'[0]==0}];
numsyseq=syseq/.{G->1,m1->2,m2->1,m3->1,r0->1}
sol=NDSolve[numsyseq,vars,{t,0,30.0}]
ParametricPlot3D[{(r[[1]]/.sol[[1]]),(r[[2]]/.sol[[1]]),(r[[3]]/.sol[[1]])},{t,0.0,30.0},PlotRange->All]