This construction looks like it should answer the question raised at the end of the previous post, but instead of checking whether it really works, I decided to post some pictures.

Given distinct points in a Euclidean space , consider the ODE system

with the initial condition . The solution exists and is well-behaved until the first collision: for some . We can stop the process then (thus obtaining a set of at most points) or remove the duplicates and restart. Let’s try the latter approach. To simplify things (at a cost), we can write the system as follows:

Now the solution exists for all times, but becomes stationary when all points come together. The cost to which I alluded above is that keeping multiple copies of collided points distorts the post-collision behavior, but I hope that the qualitative picture survives.

I used Scilab to plot the trajectories of this system for several initial configurations. For purposes of numeric computation, the ODE was modified once again:

where is tolerance, related to the step size used by the solver (Euler’s method works fine).

In each case the particles obviously want to collide: the points placed nearby move toward each other. Indeed, suppose that and are much close to each other than to the other points. The ODE system yields

where the smallness is due to the near cancellation of and for . As long as the sum of small terms is bounded by a constant less than , the distance will decrease at linear rate, heading toward zero.

Here is the (hastily written) scilab function. It takes three parameters: vector of x-coordinates, vector of y-coordinates, and the number of steps. I use 100 for the latter, increasing it if the points do not have enough time to collide. For example, try `gravity([0,1,2,3],[2,3,0,0.5],100)`

.

```
function gravity(vx,vy,Nsteps)
h = 0.01;
Npoints = length(vx);
tolerance = h*(Npoints-1)/2;
x = zeros(Nsteps,Npoints);
y = zeros(Nsteps,Npoints);
x(1,:) = vx;
y(1,:) = vy;
for m = 1:(Nsteps-1)
x(m+1,:)= x(m,:);
y(m+1,:)= y(m,:);
for i = 1:Npoints
for j = 1:Npoints
v = [x(m,j),y(m,j)]-[x(m,i),y(m,i)];
if norm(v) > tolerance then
x(m+1,i) = x(m+1,i) + h*v(1)/norm(v);
y(m+1,i) = y(m+1,i) + h*v(2)/norm(v);
end
end
end
end
clf();
plot(x(:,:),y(:,:));
endfunction
```