The N-body problem is the mathematical/dynamical problem of describing the course of multiple celestial bodies interacting via gravity (or more generally, via any type of force interacting between the bodies). The term is generally used for the "difficult" case of N greater than two, since the two-body problem has its own specific analytical solution. Three or more bodies generally requires the use of a Taylor series or N-body simulation, i.e., calculating the motions in small steps.
The case of interacting orbits has some specialized solutions: some generalized information can be calculated such as whether the orbits are likely to remain stable (a planet thrown out of the system would imply instability, for example) or become chaotic, i.e., perturbation theory.
The Jacobi integral is another method addressing three-body problems.
The million-body problem is currently (2018) barely within reach of (pure) n-body simulations but remains impractical. A specific reason for interest is that a million is the order-of-magnitude of sizeable globular clusters. Ways have been devised to make it more practical by incorporating various Monte Carlo methods, one approach being to use the MC method for weak interactions (i.e., that only slight perturbations in velocity) but simulate the strong interactions in detail. Galaxies such as the Milky Way have far more stars and are not close to being simulated in that manner. Open clusters are generally smaller, on the order of 10,000 stars, and can be simulated.