# Problem With NDsolve trying to simulate n-body-gravity problem

Im Trying to model a cloud of point masses that act according to gravity, what im strugeling with is the exclusion of cases where euclidian distance = 0, If I try to do it with an If statement in my forces function, NDsolve complains about too many dependant variables I think, but I think its rather a fault in my usage/understanding of NDsolve

Any help is appreciated

edit: im aware that the equations arent yet in working order

n = 20;
g = 1;
pos = RandomReal[{-10, 10}, {n, 3}];
masses = RandomReal[{1, 50}, n];
Tmax = 20;

Forces[xp_, yp_, zp_] :=
Sum[
masses[[i]]/
EuclideanDistance[{xp, yp, zp}, {x[i][t], y[i][t], z[i][t]}]^2
, {i, 1, n}]

Graphics3D[
Table[
Sphere[pos[[i]], ((masses[[i]]*4)/(3*Pi))^(1/3)/5]
, {i, 1, n}]
, PlotRange -> 11.5]

sol = NDSolve[
Flatten@Table[
{
(x[i])''[t] == Forces[x[i][t], y[i][t], z[i][t]][[1]]/masses[[i]],
(y[i])''[t] == Forces[x[i][t], y[i][t], z[i][t]][[2]]/masses[[i]],
(z[i])''[t] == Forces[x[i][t], y[i][t], z[i][t]][[3]]/masses[[i]],
(x[i])'[0] == 0,
(y[i])'[0] == 0,
(z[i])'[0] == 0,
x[i][0] == pos[[i, 1]],
y[i][0] == pos[[i, 2]],
z[i][0] == pos[[i, 3]]
}
, {i, 1, n}]
,
Flatten@Table[
{
x[i][t], y[i][t], z[i][t]
}
, {i, 1, n}]
, {t, 0, Tmax}]

-
You might be interested in this notebook, which contains code to deal with a fairly large number of objects specified in {mass,pos,vel} form. When two objects collide, however, a singularity arises and you simply cannot expect to integrate past that point. You can use InterpolatingFunctionDomain to deal with it gracefully but that's about it. –  Mark McClure Jun 23 '14 at 0:50