What is puzzling about MOND and to lesser extent about Dark Matter is that both uses violation of Newton mechanics to justify themselves when we know that Newtonian mechanics is wrong. The problem is that we cannot yet model behavior of galaxies using just General Relativity. Typically models assumes Newtonian physics with minimal accounting for relativity. However there is no sound proof that it is a valid approxima…
There is something that always puzzles me and I would like if something could illuminate me. How are those kind of massive scale astrophysical calculations made? Do they take into account every visible mass aggregate in the galaxy and simulate them, or just apply cinematic at scale assuming the galaxy centre is a point of certain mass? I have played a little with N body simulations and I can't but be amazed about how…
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.