I'm not much of an expert on the topic, but have also played with writing N-body simulations. At some point, obviously N^2 is unfeasible, so most simulations seem to use some approximation, such as the Barnes-Hut approach [1]. That approach divides the simulation space into an octtree. One can then use said octtree to approximate the force of distant quadrants by simply using the center of mass of an entire region. Due to the rate at which gravity drops off over large distances, it's "accurate enough", and fortunately the majority of "N" is "distant".
A good (but pretty shallow) overview exists on Wikipedia about other approaches [2] (under "Calculation optimizations").
[1] https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation
[2] https://en.wikipedia.org/wiki/N-body_simulation
EDIT: It's occurred to me that you were referring more to stability of the simulation, rather than what optimisations are used. I have less information to offer on that part, and I agree that it would be a huge problem, if complete accuracy was desired. I suspect, however, that in the kinds of simulations we're talking about, it's more about a general result such as "the whole galaxy should be divergent", not whether a single star went left or right when interacting with a black hole.