2 of 2 Corrected a math mistake, added some more explanations inspired by comments.

Numerical analysis is used to calculate approximations to things: the value of a function at a certain point, where a root of an equation is, or the solutions to a set of differential equations. It is a huge and important topic since in practice most real problems in mathematics, science and technology will not have an explicit closed-form solution (and even if they have, it might not be possible to compute with infinite precision - computers after all represent numbers with finite precision). In general there are trade-offs between accuracy and computational speed.

For the three-body problem we have three point masses in starting locations $$\mathbf{x}_i(0)$$ with velocities $$\mathbf{v}_i(0)$$ that we want to calculate for later times $$t$$. Mathematically we want to find the solution to the system $$\mathbf{x}'_i(t)=\mathbf{v}_i(t),$$ $$\mathbf{v}'_i(t)=\mathbf{f}_i(t)/m_i,$$ $$\mathbf{f}_i(t)=Gm_i \sum_{j\neq i} \frac{m_j(\mathbf{x}_j-\mathbf{x}_i)}{||\mathbf{x}_i-\mathbf{x}_j||^3}.$$

The obvious method is to think "if we move forward a tiny step $$h$$ in time, we can approximate everything to be linear", so we make a formula where we calculate the state at time $$t+h$$ from the state at time $$t$$ (and so on for $$t+2h$$ and onwards): $$\mathbf{x}_i(t+h)=\mathbf{x}_i(t)+h\mathbf{v}_i(t),$$ $$\mathbf{v}_i(t+h)=\mathbf{v}_i(t)+h\mathbf{f}_i(t).$$ This is called Euler's method. It is simple but tends to be inaccurate; the error per step is $$\approx O(h^2)$$ and they tend to build up. If you try it for a two body problem it will make the orbiting masses perform a precessing rosette orbit because of the error build-up, especially when they get close to each other.

There is a menagerie of methods for solving ODEs numerically. One can use higher order methods that sample the functions in more points and hence approximate them better. There are implicit methods that instead of trying to find a state at a later time only based on the current state look for a self-consistent late and intermediate state. Most serious methods for solving ODEs will also reduce the step-size $$h$$ when the forces grow big during close encounters to ensure that accuracy remains acceptable. As I said, this is a big topic.

However, for mechanical simulations you may in particular want to look at methods designed to preserve energy and other conserved quantities (symplectic methods - these are the ones used by professionals for long-run orbit calculations). Perhaps the simplest is the semi-implicit Euler method. There is also the Verlet method and leapfrog integration. I like the semi-implicit Euler method because it is super-simple (but being a first order-method it is still not terribly accurate): $$\mathbf{v}_i(t+h)=\mathbf{v}_i(t)+h\mathbf{f}_i(t),$$ $$\mathbf{x}_i(t+h)=\mathbf{x}_i(t)+h\mathbf{v}_i(t+h).$$ Do you see the difference? You calculate the updated velocity first, and then use it to update the positions - a tiny trick, but suddenly 2-body orbits are well behaved.

The three body problem is chaotic in a true mathematical sense. We know there are situations where tiny differences in initial conditions get scaled up to arbitrarily large differences in later positions (even if we rule out super-close passes between masses). So even with an arbitrarily fine numerical precision there will be a point in time where our calculated orbits will be totally wrong. The overall "style" of trajectory may still be correct, which is why it is OK to play around with semi-implicit Euler as long as one is not planning any space mission based on the results.