G=1 (* Stała grawitacyjna *)

ε=0.1 (* Zmiękczanie siły newtonowskiej *)

dat=Import["https://th.if.uj.edu.pl/~odrzywolek/homepage/students/aa/Plummer_16.dat","Table"]  (* 4,5,12,16,128,16384 *)

num=dat[[1,1]]

initPos=Table[dat[[ii,2;;4]],{ii,2,num+1}]

initVel=Table[dat[[ii,5;;7]],{ii,2,num+1}]

m=Table[dat[[ii,1]],{ii,2,num+1}]

(*m=Table[ToExpression["m"<>ToString[k]],{k,1,num}]*)

r=Table[{ToExpression["x"<>ToString[k]<>"[t]"],ToExpression["y"<>ToString[k]<>"[t]"],ToExpression["z"<>ToString[k]<>"[t]"]},{k,1,num}];

sys=Table[
Table[
m[[i]] D[r[[i,k]],{t,2}]==-Sum[
If[i≠j,G m[[i]] m[[j]](r[[i,k]]-r[[j,k]])/Sum[(r[[i,k]]-r[[j,k]])^2+1/3ε^2,{k,1,3}]^(3/2),0],{j,1,num}],
{i,1,num}],{k,1,3}]//Flatten;

icsPos=Table[r[[i,k]]==initPos[[i,k]],{i,1,num},{k,1,3}]/.t->0

icsVel=Table[D[r[[i,k]],t]==initVel[[i,k]],{i,1,num},{k,1,3}]/.t->0

T=25

sol=NDSolve[Join[sys,icsPos,icsVel],Flatten@r,{t,0,T}][[1]];

R=r/.sol;

ParametricPlot3D[R,{t,0,T},PlotRange->3.5{{-1,1},{-1,1},{-1,1}},PlotStyle->Thin]

Manipulate[ListPointPlot3D[R/.t->τ,PlotRange->4{{-1,1},{-1,1},{-1,1}},PlotStyle->Thin,BoxRatios->{1, 1, 1}],{τ,0,T,0.1}]