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=2;
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};
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,Table[(r[[1,j]]/.{t->0})==0,{j,1,3}],Table[(D[r[[1,j]],{t,1}]/.{t->0})==0,{j,1,3}],{x2[0]==1.0,y2[0]==1.0,z2[0]==0,x2'[0]==-10.0,y2'[0]==0,z2'[0]==0}]

numsyseq=syseq/.{G->1,m1->1000,m2->1,m3->1}
sol=NDSolve[numsyseq,vars,{t,0,0.5}]
ParametricPlot3D[{(r[[1]]/.sol[[1]]),(r[[2]]/.sol[[1]])},{t,0.0,0.5},PlotRange->{{-0.2,1.2},{-0.2,1.2},{-0.1,0.1}}]

vals={G->100,m1->10,m2->5,m3->1};
numsyseq=syseq/.vals;
sol=NDSolve[numsyseq,vars,{t,0,0.5}]
ParametricPlot3D[{(r[[1]]/.sol[[1]]),(r[[2]]/.sol[[1]])},{t,0.0,0.5},PlotRange->All]

rcm=(m1 r1+m2 r2)/(m1+m2)/.vals/.sol[[1]];
r1r=r1-rcm/.sol[[1]];
r2r=r2-rcm/.sol[[1]];
ParametricPlot3D[{r1r,r2r},{t,0,0.5}]