The main issue - and difference - between the $N$-body simulations that we are capable of running on computers and the gravitational interactions in the real universe is that we discretise time into a series of time-steps in order to "step" the simulation forward.
On the other hand, in the universe, time is continuous, and is not "chopped" up into intervals of any kind. The movement of bodies is continuous through space. Therefore, when two or more bodies are interacting under gravity somewhere in the cosmos, the force - and, therefore, acceleration - experienced by each body at every moment in time changes continuously with their changing position.
In our $N$-body simulations, however, since time is discretised into a series of steps, the movement of bodies is approximated by only evaluating the forces on bodies at the beginnings and ends of time-steps. The consequence of this is an error in the precise positions (and velocities) of the bodies, since the force on the bodies was not evaluated at the infinitely many points between the times $t_n$ and $t_{n+1}$, when the force was actually changing (and so the effect of these changing forces was not reflected at all in the positions and velocities of the bodies between $t_n$ and $t_{n+1}$).
Of course, there is little that can be done about this, since, although the size of the time-step, $\Delta{t}$, can be made arbitrarily small (albeit at the cost of computation speed), it will still be infinitely larger than zero, and thus errors will always begin to accumulate.
On top of this, as you have alluded to, there is no closed-form, analytic solution to the general $N$-body problem (for $N>2$) - hence why we resort to numerical methods - so, for now, we are stuck with our current simulations! That being said, a lot of them are extremely accurate and/or impressive in the amount of bodies they can evolve based on a variety of different integration methods and/or force approximation techniques.