for (let i = 0; i < stars.length; i++) {
    for (let j = 0; j < stars.length; j++) {
      // We don't want to calculate the same star against itself
      if (i === j) continue;

Notice this processes each pair twice: First i,j and later j,i again, doing the same calculations.

You can just write:

for (let i = 0; i < stars.length-1; i++) {
    for (let j = i+1; j < stars.length; j++) {

Does only half of the work, and no self-check needed.

But then you also need to change the Attract function so it adds the force to both bodies. Thus it should be no member function of the star class, as it affects a pair of stars.
In practice it should run twice as fast after the changes.

jonny-b said:
There is a black hole star that has more mass than the others combined and it doesn't interact in the way I thought.

It has just 10 times the mass, that's not much. I tried 1000 times. The whole universe started to wobble in a interesting way.
But the distance clamping will surely have an effect. I do not see a need to constrain max distance which causes only weak forces. You could just constrain the minimum distance to avoid div by 0.
Or add some constant (to prevent discontinuities from clamping), which i have tried:

// dsquared = constrain(dsquared, 100, 9999);
    dsquared += 100;

I also put the big star to the center, so i can more likely see what happens to that:

// Black hole
  stars.push(new Star(750/2, 750/2, 100000, 1));

But it flies out of the canvas, whil i would expect it sticks close to the center, attracting all other stars.

Maybe you did not respect star mass when integrating the force…

…yes:

update() {
    this.vel.add(this.acc); // it is force, not yet acceleration. To convert it to acceleration we need mass: acc = force / mass
    this.pos.add(this.vel);
    this.acc.mult(0);
  }

I fix this:

update() {
    this.vel.add(this.acc / this.mass);
    this.pos.add(this.vel);
    this.acc.mult(0);
  }

Now math may be right, but variable name is wrong. You need to change it to force, or do the division in the Attract function already, but that's n*n divisions vs. n divisions, so slower.

However, now stuff does not move at all, or very slowly.
Tried to increase masses but no luck. Another option would be to introduce a timestep larger than one. You can see this in my code example above.

jonny-b said:
Both of them are exerting a force on the other and we need to capture that.

Yes, but forceAtoB = -forceBtoA.

jonny-b said:
In your model I think one of those stars doesn't get tugged on.

Make a little script and print both i and j from the inner loop. You will see it covers all possible pairs, skips i==j, and does no pair twice. Do this to eliminate your doubt.

It's no optimization, just the least you should do for n-body problems. If you don't, the impression is you can not write performant code. /:O\

Otherwise the simulation looks nice now : )

jonny-b said:
Hmm, I did your solution with just a black hole and a single star and the black hole never moved at all. Then I tried all stars and I noticed that one star would never move.

Then you did something wrong. It's not ‘my’ method. Everybody uses this, there are just some variants of this same thing.

A matter of minutes to proof it, but do i need to register to fork your project, and be able to save my changes?

No problem, but you must be stopped doing this ;D

Seems ok: https://codepen.io/joedoej/pen/xxYPjJZ?editors=0011​

Now you can do twice the stars in the same time.