num=3

m={1,2,4}

r=Table[{ToExpression["x" <> ToString[i]<>"[τ]"],ToExpression["y" <> ToString[i]<>"[τ]"],ToExpression["z" <> ToString[i]<>"[τ]"]},{i,1,num}]

f=Table[Sum[If[ii==jj,
0,
-m[[ii]] m[[jj]](  r[[ii]]-r[[jj]]  )/(Sqrt[(r[[ii]]-r[[jj]] ).(r[[ii]]-r[[jj]])]^3)],{jj,1,num}],{ii,1,num}];

f[[2,3]] ;

syss=Table[{m[[i]]D[ToExpression["x" <> ToString[i]<>"[τ]"],τ,τ]==f[[i,1]],m[[i]] D[ToExpression["y" <> ToString[i]<>"[τ]"],τ,τ]==f[[i,2]],m[[i]] D[ToExpression["z" <> ToString[i]<>"[τ]"],τ,τ]==f[[i,3]]},{i,1,num}]//Flatten;

sys=syss;

funs=Table[{ToExpression["x" <> ToString[i]<>"[τ]"],ToExpression["y" <> ToString[i]<>"[τ]"],ToExpression["z" <> ToString[i]<>"[τ]"]},{i,1,num}]//Flatten

ics=Table[{ToExpression["x" <> ToString[i]<>"[0]"]==RandomReal[{0,1}],ToExpression["x" <> ToString[i]<>"'"<>"[0]"]==RandomReal[{-1,1}],
ToExpression["y" <> ToString[i]<>"[0]"]==RandomReal[{0,1}],
ToExpression["y" <> ToString[i]<>"'"<>"[0]"]==RandomReal[{-1,1}],
ToExpression["z" <> ToString[i]<>"[0]"]==RandomReal[{0,1}],
ToExpression["z" <> ToString[i]<>"'"<>"[0]"]==RandomReal[{-1,1}]},{i,1,num}] //Flatten

problem=Union[{sys,ics}];

sol=NDSolve[problem,funs,{τ,-3.3,3.3},MaxSteps->200000]

traj=r /. sol //First

ParametricPlot3D[{traj[[1]]-(m[[1]]traj[[1]]+m[[2]]traj[[2]]+m[[3]]traj[[3]])/(Tr[m]),traj[[2]]-(m[[1]]traj[[1]]+m[[2]]traj[[2]]+m[[3]]traj[[3]])/(Tr[m]),traj[[3]]-(m[[1]]traj[[1]]+m[[2]]traj[[2]]+m[[3]]traj[[3]])/(Tr[m])},{τ,-3.3,3.3
},PlotRange->{{-0.91,0.91},{-0.91,0.91},{-0.91,0.91}},PlotStyle->{Red,Blue,Black},AspectRatio->1,PlotPoints->1000]

t1=Table[traj[[1]]-(m[[1]]traj[[1]]+m[[2]]traj[[2]]+m[[3]]traj[[3]])/(Tr[m]),{τ,-3.3,3.3,0.001}];

t2=Table[traj[[2]]-(m[[1]]traj[[1]]+m[[2]]traj[[2]]+m[[3]]traj[[3]])/(Tr[m]),{τ,-3.3,3.3,0.001}];

t3=Table[traj[[3]]-(m[[1]]traj[[1]]+m[[2]]traj[[2]]+m[[3]]traj[[3]])/(Tr[m]),{τ,-3.3,3.3,0.001}];

Animate[ListPointPlot3D[{{t1[[ii]]},{t2[[ii]]},{t3[[ii]]}},PlotRange->{{-1,1},{-1,1},{-1,1}},PlotStyle->{Red,Blue,Black},AspectRatio->1],{ii,1,Dimensions[t1][[1]],1}]