N Gravitational Body Simulation I'm trying to compute an N-Body simulation in Matlab, where I have randomly generated 100 points in 2 dimensions and directly computed the acceleration formula. I believe by doing this I have done it in specifically one time frame, and I want to extend it to multiple time frames but I'm not really sure how...? I have posted the code below and any help of where to go next would appreciated. (Note, in this case I have assumed that the mass M is constant (that of the milky way). The end result it to key in the 127 nearest galaxies to the milky way to see how the galaxies interact but I wanted to see how to code a simple case first)
My code doesn't seem to follow an ODE method like I have seen elsewhere so if my code is completely wrong I'd also be grateful to know (I'm not very experienced in Matlab)!! Thanks.
clear all
N=100; %Number of galaxies
rangex = 100;
rangey = 100;
T=100; %Number of timesteps
G = 6.67408*(10^-11); %Gravitational constant in m^3 kg^-1 s^-2
M = 1.15*1.98855*(10^42); %Mass of MW in kg

for t=1:T
for i=1:N
    x(i) = (rand-0.5)*rangex;
    y(i) = (rand-0.5)*rangex;
    P{i}=[x(i),y(i)];
end

for i=1:N    

    xs(i) = 0;

    for j=1:i-1
        dx(i) = P{j}(1,1)-P{i}(1,1);
        xs(i)=xs(i)+G*M*((P{j}(1,1)-P{i}(1,1))/(sqrt((P{j}(1,1)-P{i}(1,1)).^2)).^3);
    end

    for j=i+1:N
        dx(i) = P{j}(1,1)-P{i}(1,1);
        xs(i)=xs(i)+G*M*((P{j}(1,1)-P{i}(1,1))/(sqrt((P{j}(1,1)-P{i}(1,1)).^2)).^3);
    end

end
for i=1:N    

    ys(i) = 0;

    for j=1:i-1
        ys(i)=ys(i)+G*M*((P{j}(1,2)-P{i}(1,2))/(sqrt((P{j}(1,2)-P{i}(1,2)).^2)).^3);
    end

    for j=i+1:N
        ys(i)=ys(i)+G*M*((P{j}(1,2)-P{i}(1,2))/(sqrt((P{j}(1,2)-P{i}(1,2)).^2)).^3);
    end
end

end

quiver(x,y,xs,ys)

 A: Consider the formula you were given, which (based on a comment)
was equivalent to this:
$$
−G \sum_{\substack{1\leq j\leq N \\[.2ex] j\neq i}}
  m_j \frac{(\vec x_j − \vec x_i)}{(\lVert \vec x_j − \vec x_i\rVert)^3 }
$$
Where you had $\mathrm{abs}(\cdots)$,
I've written $\lVert \cdots \rVert$ because that's a typical way of
denoting the magnitude of a vector.
I also noticed that the iteration should be over $j$ for a constant
$i$ (so that $m_j$ has to be inside the summation, unlike the
constant factor of $G$, which can be moved outside);
but I think you got that part right in your implementation.
This formula is essentially one side of the differential equation
of motion of the galaxies in your simulation. 
That's usually all you have in a simulation like this. 
But the little arrows on top of variables such as $\vec x_i$ means they are vectors, that is, $\vec x_i$ represents both coordinates of one of the galaxies, not just the "$x$" coordinate.
The numerator inside the summation is a simple coordinate-wise sum,
so in principle you could compute it for all the coordinates on one
axis first, and then for all the coordinates on the other axis.
But the denominator cannot be done that way.
Suppose you had two galaxies whose coordinates differed by a billion light-years in one dimension, but by only one million light-years in the other dimension.
The total force of one galaxy on the other is determined by the
billion-light-year distance, but when you look at the second coordinate,
you're only using one million light-years as the separation, so you're computing a force that is at least a million ($1000^2$) times as powerful as it should be.
Rather than writing the denominator in the
form $(\lVert \vec x_j − \vec x_i\rVert)^3$, try
$$
\left(\vec r^T \; \vec r\right) ^{1.5},
$$
where $\vec r = \vec x_j − \vec x_i.$
The reason this works is that the inner product (aka dot product)
of a vector with itself is the square of the vector's magnitude;
and as a matrix operation, the inner product of a column vector
is written as the product of the transpose of the vector on the left and the vector itself on the right. That is,
$\vec r^T \; \vec r = (\lVert \vec x_j − \vec x_i\rVert)^2$. 
Rather than take the square root of this and then cube it, we can equivalently just take the $3/2$ power.
In vector notation, then, the formula is
$$
−G \sum_{\substack{1\leq j\leq N \\[.2ex] j\neq i}}
  m_j \frac{\vec r}{\left(\vec r^T \; \vec r\right) ^{1.5}},
\quad \text{where } \vec r = \vec x_j − \vec x_i.
$$
MATLAB gives really good support for matrix and vector operations,
so you might as well use them rather than trying to pick apart the
components of the vectors.
As for starting with a simple example, a really good place to start
would be with just two galaxies. Put them at some distance from each other
and give each one a starting velocity so that they will both orbit in the
same direction around their common center of mass.
The two-body problem is so thoroughly solved that you should be able to
find the orbital speed that will make the two orbits circular.
So set up those initial conditions, let the simulation run for a
few orbits, and see if the bodies maintain circular orbits of
constant radius, as they should.
(Personally, it's been decades since I tried to do this, but I never
had much luck with the obvious Euler-method simulation,
which is why I suggest to start with a scenario with a
known solution like this before tackling the harder stuff.)
