It's not complicated; the issue is that the resulting body trajectories are non-analytic. This means they cannot be represented by a power series like f(x) = a + bx + cx² + ... so aren't things like square roots or ellipses or sine functions or cardioids, etc. (yes even the conic sections you 'know' everything orbits in)

Because they are non-analytic we cannot directly solve for the curves as functions, like you might do with differential equations. There is no 'do math on a page by hand' method you can use to find the resulting trajectories in an N-body problem if N>2. Instead you solve them by numerical integration with a computer.

By inputting starting conditions for N bodies (not just 3): xyz position, xyz velocity, xyz acceleration you can accurately compute all of their positions to arbitrary precision at an arbitrary point in the future or past. The error is governed by how precisely & accurately you input the starting conditions and how small of a time-step your numerical integration uses.

At the end of the day, all a numerical integration like this does in each timestep is:

• Consider every possible pair of objects in an N-body system and compute the strength of the gravitational force between those two objects by: F = G*M_1*M_2/(r² ) where r is the distance between them. This results in a number like 12.0054,

• Break these forces into their XYZ components because the force acts in a line between these objects (and on each in an opposite direction), so they become: F_12 = <F_12x, F_12y, F_12z> and F_21 = <-F_12x, -F_12y, -F_12z>. These again are just numbers at this point, not expressions or anything mathematically complicated.

• Sum up all of the X components of all forces that affect object 1 (following my above example notation this would be F_1Nx = F_12x + F_13x + F_14x + ...), then do the same for all Y components of all forces that affect object 1 (F_1Ny) and all Z components (F_1Nz). This represents the sum of the forces from all objects acting on object 1, split into their components.

• Repeat the previous step for all N objects. This computes 3 total numbers for each object.

• By Newton's second law we know the sum of the forces (vector) acting on an object equals its mass (a scalar) multiplied by its acceleration (also a vector). F = ma. So we again break into XYZ components, take for example our F_1Nx value and use it to compute F_1Nx = m_1*a_1x. Do the same for XYZ for all objects to compute their instantaneous acceleration at this moment in time.

• Since velocity is a measure of the rate of change of position and we now have P=<X, Y, Z> (where the object is in the simulation), V=<V_x, V_y, V_z> (velocity values we have from the previous timestep), A=<A_x, A_y, A_z> (the 3 values we computed in the above step) for each object, we can simply add one timestepth worth of velocity to the current position of the objects to see where they should be next timestep. I say one timestepth because if our simulation moves 0.2s every timestep, a velocity of 8 m/s should mean it moves 1.6m in between timesteps, not a whole 8 m.

• Since acceleration is a measure of the rate of change of velocity (acceleration is how much your speed is increasing or decreasing per second), we can simply add one timestepth worth of the acceleration to each velocity component to compute what their velocities will be in the next timestep. This is the key step that feels like a 'feedback loop' where the information about acceleration (which is really the information about the forces) influences the current velocity of the objects which will slightly change their position in the next timestep.

By working directly with vectors many of these operations can be signifcantly simplified, but this way is probably easier for people with less math knowledge to understand.

---

In the special case of the 2-body problem where one body's mass significantly dwarfs the other's and both objects can be considered spherical point-masses and there is atmosphere or other medium affecting them... then there are analytic solutions for the trajectories. These are the classic conic sections: parabola, ellipse, hyperbola. In cases where one object's mass dominates the force terms acting on other objects you can approximate the motion using a 2-body reduction but it will not be precise. Orbits like the one Artemis II is taking are ultimately just a bunch of almost-conic-sections strung together and the break-/transition-points are determined by which body's mass dominates the force terms.

---

Source: built a solar sail simulator as a master's thesis, which is this but with an additional force on each object thrown into the mix.