Formally, if gr in some situations approaches or becomes newtonian theory of gravity, then necessarily there is no known analytic solution to the N body problem. For 2 body the eccentric anomaly (typically a variable ultimately sought when trying to find space and time parameterization of the motion of orbiting bodies) can be solved for exactly though usually only after an initial guess and the methods differ for different eccentricities. Ultimately what happens is you are trying to solve an equation of the form M=E+asinE where M is the known mean anomoly and E the eccentric anomaly, for E.
This is considered by many phisicist from Kepler, Newton to current day an important and difficult question. It is where many of the fields of chaos thoery have come from. Basically the study of systems where very small pertubations to the initial conditions give drastically different outcomes. This happens with multiple body's all having central force fields (ie gravity, electric fields). Unfortunately the problem does not go away with General Relativity but wrather increases in difficulty.
Frankly a Relativistic N body simulation that is accurate to a certain degree (converges to actual solution closely enough) would be very difficult to write and is probably considered "state of the art" for both modern and classical orbital mechanics and computer science, as it is computationally difficult for computers.
One way to attempt a solution would be to compute all the relativistic dynamics of the initial conditions and use the principle of superposition, then neglecting the gravitational fields, compute a small time interval "step". Then you would continue doing this. There would be significant flaws however as I said earlier n body problem is chaotic and slight error will be magnified and only a few iterations in the model could be completely unrealistic. Also tensors are a must for implementing and understanding general relativity, also tensors (matrices) will ease computation and code writing, hence their development.
Also guages are a topic of partial differential equations, with out going too far into them (you really would have to read a text on pde's), they are similar (i dare say) to the concept of a constant of integration. basically a constant or variable which was not in the original problem but turns up when solving the pde's. they are basically a vector, matrix or tensor though not a scalar and interact with differential operators (also the reason for them existing exactly what they are). They are used to solve pde's as you can sometimes make them equal anything and basically change the pde so that you can solve it.
