Below is a code, in python, that gives the position of two particles in a circular orbit. All constants are set to 1 for now to make it simple, and will be changed when real data is to be plugged in. The code needs to be modified to allow a third particle in the system and give the positions of the third particle as well as the first two. Any notes as to the changes made would be much appreciated so it can be more easily understood how it works.

import numpy as np from scipy.integrate import ode G = 1 n = 2 def initialize_binary_circular_orbit(): m = np.zeros(n) y = np.zeros(6 * n) m[0] = 1 #mass of particle 1  m[1] = 1 #mass of particle 2  a = 1 #distance between particle 1 and 2  omega = np.sqrt(G * (m[0] + m[1]) / a ** 3) #Newton's Third Law "Angular Velocity"  y[0] = a * m[1] / (m[0] + m[1]) # x for particle 1  y[1] = 0 # y for particle 1  y[2] = 0 # z for particle 1  y[3] = 0 # vx for particle 1  y[4] = omega * y[0] # vy for particle 1  y[5] = 0 # vz for particle 1  y[0 + 6] = -a * m[0] / (m[0] + m[1]) # x for particle 2  y[1 + 6] = 0 # y for particle 2  y[2 + 6] = 0 # z for particle 2  y[3 + 6] = 0 # vx for particle 2  y[4 + 6] = omega * y[0 + 6] # vy for particle 2  y[5 + 6] = 0 # vz for particle 2  return m, y def f(t, y): dydt = np.zeros(6 * n) # d(position)/d(time)  for i in range(n): vx = y[3 + 6 * i] vy = y[4 + 6 * i] vz = y[5 + 6 * i] dydt[0 + 6 * i] = vx dydt[1 + 6 * i] = vy dydt[2 + 6 * i] = vz # d(velocity)/d(time)  for i in range(n): xi = y[0 + 6 * i] yi = y[1 + 6 * i] zi = y[2 + 6 * i] for j in range(n): if i >= j: continue  xj = y[0 + 6 * j] yj = y[1 + 6 * j] zj = y[2 + 6 * j] r = np.sqrt((xi - xj) ** 2 + (yi - yj) ** 2 + (zi - zj) ** 2) Fx = G * m[i] * m[j] * (xi - xj) / r ** 3 # x-component of force on j due to i  Fy = G * m[i] * m[j] * (yi - yj) / r ** 3 # y-component of force on j due to i  Fz = G * m[i] * m[j] * (zi - zj) / r ** 3 # z-component of force on j due to i  dydt[3 + 6 * i] = dydt[3 + 6 * i] - Fx / m[i] dydt[4 + 6 * i] = dydt[4 + 6 * i] - Fy / m[i] dydt[5 + 6 * i] = dydt[5 + 6 * i] - Fz / m[i] dydt[3 + 6 * j] = dydt[3 + 6 * j] + Fx / m[j] dydt[4 + 6 * j] = dydt[4 + 6 * j] + Fy / m[j] dydt[5 + 6 * j] = dydt[5 + 6 * j] + Fz / m[j] return dydt # main program here  m, ystart = initialize_binary_circular_orbit() tstart = 0.0 r = ode(f, jac=None).set_integrator('dopri5') r.set_initial_value(ystart, tstart) tend = 1.e1 dt = 0.1 while r.successful() and r.t < tend: r.integrate(r.t + dt) print (r.t,r.y[0],r.y[1],r.y[6],r.y[7]) # r.t is tfinal, r.y[0] is x1, r.y[1] is y1, r.y[6] is x2, r.y[7] is y2