The key to simulating the motions of bodies in solar system is the `scipy.integrate.odeint()` method. This function takes three arguments: (a) an initial state, (b) an array of variable values to evaluate the system at (in this case, these are instants of time) and (c) a function `ode()` that, given the values that characterize the state at one point returns the derivatives for those values.

Since we need the derivatives for all the values that characterize the solar system, initial_state encodes the state of the solar system in one long list of 16*7 = 112 items. (16 objects, 3 position and 3 velocity vectors + mass for each object.) The `planet_seperator()` function seperates this long list into individual lists for each body.

How do we write `ode()`? Well,the derivatives for the position vectors (x,y,z) are just the velocity vectors (vx, vy, vz). The acceleration vectors are calculated according to Newton's Law of universal gravitation. For each pair of bodies, o1 and o2, the component acceleration of o1 due to o2 is G (the gravitational constant) times the mass of o2 divided by the square of the distance between o1 and o2. The resultant acceleration is the sum of all the component accelerations.