[Numerical analysis][1] 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][2]. 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)=\sum_{j\neq i} Gm_i m_j/||\mathbf{x}_i-\mathbf{x}_j||.$$

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][3]. 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][4] 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][5] 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. As I said, this is a big topic.

However, for mechanical simulations you may want to look at methods designed to preserve energy and other conserved quantities ([symplectic methods][6] - these are the ones used by professionals for long-run orbit calculations). Perhaps the simplest is the [semi-implicit Euler method][7]. There is also the [Verlet method][8] and [leapfrog integration][9]. 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.

  [1]: https://en.wikipedia.org/wiki/Numerical_analysis
  [2]: https://en.wikipedia.org/wiki/Closed-form_expression
  [3]: https://en.wikipedia.org/wiki/Euler_method
  [4]: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
  [5]: https://en.wikipedia.org/wiki/Explicit_and_implicit_methods
  [6]: https://en.wikipedia.org/wiki/Symplectic_integrator
  [7]: https://en.wikipedia.org/wiki/Semi-implicit_Euler_method
  [8]: https://en.wikipedia.org/wiki/Verlet_integration
  [9]: https://en.wikipedia.org/wiki/Leapfrog_integration