Astronomers have almost half a century of experience in writing computer codes to simulate dense stellar systems. The first published results date back to 1960, and it was in the subsequent decade that it became clear just how tricky it was to simulate a group of interacting stars. The task seems so easy: just integrate Newton's equations of motion for each star, under the influence of the gravitational pairwise interactions of all other stars. Indeed, it is straightforward to write a simple code to do so, as we will see below. And as long as all stars remain fairly well separated from each other, even a simple code will do a reasonably good job.
In practice, though, even a small group of stars will spontaneously form one or more double stars. This was discovered experimentally in the early sixties. One way to understand this result, after the fact, is from an energetic point of view. When a double star, or binary as they are generally called, is formed, energy is released. Since the two stars in a binary are bound, the potential energy is larger than the kinetic energy, and the total energy is negative. When three stars come together randomly, there is a chance that two of the three are left in a bound state, while the third one escapes, carrying the excess energy. Left by itself, a stellar system will exploit this energy liberation mechanism by spontaneously forming binaries.
As soon as even one binary appears, a simple code with constant time steps will either crash or resort to a very slow and inefficient crawl. The changes in the relative configuration of the two stars occur so frequent that they can be easily missed. If they are modeled correctly, by choosing a tiny time step, the rest of the system will be forced to slow down, while spending a large and unnecessary amount of computer time on all the other stars that don't need such fine resolution. By the end of the sixties, these problems were overcome by the development of codes that employed individual time steps. Stars with close neighbors were stepped forward in time more frequently than stars at large, and in this way the computational power was brought to where it was most needed.
This modification in itself brought gravitational 
-body codes already
well outside the range of systems that are normally discussed in text
books on numerical integration methods. The internal book keeping
needed to write a correct and efficient code with individual time
steps is surprisingly large, given the simplicity of the task:
integrate the effect of pairwise attractive inverse square forces.
But this was only a first step toward the development of modern 
-body
codes. The presence of tight binaries produced much more of an obstacle,
and throughout the seventies a variety of clever mechanisms were developed
in order to deal with them efficiently.
For one thing, there are problems with round-off. Two stars in a tight orbit around each other have almost the same position vector, as seen from the center of a star cluster, where we normally anchor the global coordinate system. And yet it is the separation between the stars that determines their mutual forces. When we compute the separation by subtracting two almost identical spatial vectors, we are asking for (numerical) trouble. The solution is to introduce a local coordinate system whenever two or more stars undergo a close interaction. This does away with the round-off problem, but it introduces a host of administrative complexities, in order to make sure that any arbitrary configuration of stars is locally presented correctly - and that the right thing happens when two or more of such local coordinate patches encounter each other. This may not happen often, but one occurrence in a long run is enough to crash the system if no precautions have been taken for such a situation to happen.
We can continue the list of tricks that have been invented to allow every larger and denser systems to be modeled correctly. Numerical problems with the singularity in the two-body system have been overcome by mapping two or more interacting stars from the three-dimensional Kepler problem to a four-dimensional harmonic oscillator. The total force on particles has been split into different contributions, the first from a near zone of relatively close neighbors and the second from a far zone of all other particles, with each partial force being governed with different integration time steps. Tree codes have been used to group the contributions of a number of more and more distant zones together in ever larger chunks, for efficiency. Triple stars have received their own special treatment, especially the marginally stable triples that are sometimes long-lived, but continuously changed their inner state due to internal perturbations. The list goes on.
In this first book, we will introduce a modern integrator, the Hermite scheme, developed in the 1990s, together with a variable time step integration scheme, where all stars share a common time step at any given time. In subsequent volumes in the series, we will introduce the other refinements mentioned above. Our emphasis will be on a complete explanation of all the steps involved, together with a discussion of the motivation for those steps. In the last few chapters, we will embark on a research project featuring stellar collisions, in a simple gravity-only approximation.