In a particle simulation, it's pretty easy to compute gravitational interactions directly from Newton's Law of Universal Gravitation: f = g m1 m2 / r

g is the gravitational constant

m1 and m2 are the masses of the objects

r is the distance between the centers of the objects

f is the magnitude of the force

Since you want a vector force, you need to scale this magnitude f by a vector, which can be the normalized direction vector R between the object centers. Or you can save one vector scale factor by using the non-normalize vector between centers R, and calculating F = R * g m1 m2 / r

The other problem in a simulation is particles will occasionally pass very close to one another, resulting in a distance r near zero, which gives a force near infinity. To avoid this, the force is typically "softened" by dividing by 1/(a+r

See my simple N

One fairly simple way to get much better performance is to only consider nearby neighbors. You can find your nearby neighbors by storing lists of particles by their integer coordinates in a hash table, and looking through your list to find your closest neighbors.

See my hashtable
gravity simulation, which immediately scales to over 1000
particles. But there's a problem: the galaxy can't hold
together anymore, even if we artificially inflate the
gravitational constant, because the long-range forces that bind
one side of the galaxy to the other are never evaluated.

The hashtable version gives the wrong answer because it ignores
long-range interactions. There are several neat ways to
approximate these interactions.

A very simple approach is to sum
all the particles into one "megaparticle", which will end up
at the center of the galaxy. We can then compute each
particle's interaction with its local neighbors, and then with the
megaparticle. This doesn't even conserve momentum, because
the influence of a particle on the megaparticle isn't symmetric,
but at least the galaxy holds together!

A more accurate version of this is called Barnes-Hut:
for neighboring particles, we compute full particle-particle
interactions. To interact with a group of distant particles,
we lump the group into a megaparticle. To accurately define
which particles form a "group", we impose a spatial tree across
the simulation, using the coarse levels of the tree to approximate
distant regions, and increasingly fine levels for more nearby
regions.