After Newton discovered the equations behind Kepler's 2-Body solution, he ran afoul of the n-Body Problem. A limited solution presented itself in the equilibrium point; the full five solutions were offered by Lagrange. Further refinements handled special cases of the 4-body like when Laplace sussed out three of the heaviest (and largest) Jovian moons. So here we are in the computer age, when computers struggle to handle particles numbering more than maybe five. It's O(n*n).
Kepler and Newton assumed that the bodies are rigid points and that we ignore Einstein; so Io's tides and Mercury's depth in the well don't figure. Also computers have to approximate floating point numbers well above Zeno's quantum level (Planck?). That means we have already sold the horse for approximations. And that means we may as well approximate some more. Keyframe Codes has a little primer on the available algos.
Midlevel nerds have been using Barnes-Hut. This method runs on multipoles: it buckets the particles in squares-of-particles, finds the center of gravity, and runs the calc for the centers. The buckets are identified with quad-trees; Barnes-Hut ends up O(n * ln(n)). Leave aside for now if we want to run this 2D or 3D; we'll pretend perfect ecliptic. What if the cluster doesn't behave like a square (or cube). Real astrodynamicists need the particle swarm where most particles would be clustering. They will be looking first at Lagrange's L4 and L5 "Trojan" haloes especially for Jovians. These haloes are elongates.
In 1987 Leslie Greengard came up with something else, the "Fast Multipole Method" - usually FMM. Not only does it claim O(n) but, depending on how often you want to run each n, the error can arbitrarily diminish. This was Greengard's thesis and, it seems, got a refinement under one Vladimir Rokhlin. Discontents exist: Srinivas Aluru in "Greengard’s N-body Algorithm is not order N" (pdf) showed that the choice of precision must depend on n (at least ln n), therefore cf the title - to the tune of O(n ln2n). That extra ln(n) is, then, the price for accuracy Barnes-Hut must pay in the Trojan haloes. Dr Greengard's peers seem to have agreed; he got to keep his PhD and the IEEE considers FMM in the top ten of twentieth-century algos.
This starts by keeping each interaction as a field vector Φ and treating that vector as a number on the complex plane. The vector acting upon particle "i" is the gravitational parameter of counterpart "j", times the natural logarithm of the i-j distance (plus iθ); that's the force. The iθ exists, here, to make it a vector as to place the particle on the coördinate plane; a computer may ignore it now to figure the coörds when it's display-time.
Some might then use "numerical methods", which is mathematician for iteration. The FMM, in the complex plane, may use the Cauchy-Riemann equations. At this point I'm leaving the rest to the nerds.