Plot orbit of three bodies using rk4

I have been trying to plot the trajectories of three particles using the RK4 method. I haven't been able to produce an array of results over the time period as it brings up the following error message:

  File "C:\Users\Local\Runge-Kutta 4 Code.py", line 65, in <module>
    solution.step()

  File "F:\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 170, in step
    raise RuntimeError("Attempt to step on a failed or finished "

RuntimeError: Attempt to step on a failed or finished solver.

I suspect that there is a problem with the initial "y_0" value that I have but I could be wrong. Any help at all would be greatly appreciated. My code is as follows:

import numpy as np
import matplotlib.pyplot as plt
import scipy as scipy

from scipy import integrate
from numpy import asarray
from numpy import savetxt

# Physical constants
mass_vector = np.array([1, 1, 1])

r_vec_1 = np.array([0, 0])
v_vec_1 = np.array([-np.sqrt(2), -np.sqrt(2)])

r_vec_2 = np.array([-1, 0])
v_vec_2 = np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])

r_vec_3 = np.array([1, 0])
v_vec_3 = np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])


#Initial x acceleration ODE's
def x1_double_dot(y, mass_vector):
    #G can be omitted for scale purposes (should not be compared with realistic data)
    return ((mass_vector[1]*(y[4]-y[0])/((y[0]-y[4])**2 + (y[1]-y[5])**2)**(3/2)) +
            (mass_vector[2]*(y[8]-y[0])/((y[0]-y[8])**2 + (y[1]-y[9])**2)**(3/2)))
def x2_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[0]-y[4])/((y[4]-y[0])**2 + (y[5]-y[1])**2)**(3/2)) +
            (mass_vector[2]*(y[8]-y[4])/((y[4]-y[8])**2 + (y[5]-y[9])**2)**(3/2)))
def x3_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[0]-y[8])/((y[8]-y[0])**2 + (y[9]-y[1])**2)**(3/2)) +
            (mass_vector[1]*(y[4]-y[8])/((y[8]-y[4])**2 + (y[9]-y[5])**2)**(3/2)))
#Initial y acceleration ODE's
def y1_double_dot(y, mass_vector):
    return ((mass_vector[1]*(y[5]-y[1])/((y[0]-y[4])**2 + (y[1]-y[5]**2)**(3/2))) +
            (mass_vector[2]*(y[9]-y[1])/((y[0]-y[8])**2 + (y[1]-y[9])**2)**(3/2)))
def y2_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[1]-y[5])/((y[4]-y[0])**2 + (y[5]-y[1]**2)**(3/2))) +
            (mass_vector[2]*(y[9]-y[5])/((y[4]-y[8])**2 + (y[5]-y[9])**2)**(3/2)))
def y3_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[1]-y[9])/((y[8]-y[0])**2 + (y[9]-y[1]**2)**(3/2))) +
             (mass_vector[1]*(y[5]-y[9])/((y[8]-y[4])**2 + (y[9]-y[5])**2)**(3/2)))

#This is my X(t) at time zero
y_0 = np.concatenate((r_vec_1, v_vec_1, r_vec_2, v_vec_2, r_vec_3, v_vec_3))
y = y_0

#This is my F(X) at time zero
def fun(t,y):
    return np.array([y[2], y[3], x1_double_dot(y, mass_vector), y1_double_dot(y, mass_vector),
              y[6], y[7], x2_double_dot(y, mass_vector), y2_double_dot(y, mass_vector),
              y[10], y[11], x3_double_dot(y, mass_vector), y3_double_dot(y, mass_vector)])

# collect data
t_values = []
y_values = []

#Time start, step, and finish point
t0,tf,t_step = 0, 2, 0.1
nsteps = int((tf - t0)/t_step)

solution = integrate.RK45(fun, t0, y_0, tf, first_step=t_step)

#The loop for running the Runge-Kutta method over some time period.
for step in range(nsteps):
    solution.step()
    y_values.append(solution.y[0])
    # break loop after modeling is finished
    if solution.status == 'finished':
        break

I condensed the setup and derivatives computation to

masses = [1,1,1]
r1,v1 = [ 0,0], [-2**0.5,-2**0.5]
r2,v2 = [-1,0], [0.5**0.5, 0.5**0.5]
r3,v3 = [ 1,0], [0.5**0.5, 0.5**0.5]
G = 1

def odesys(t,u):
    def force(a): return G*a/sum(a**2)**1.5
    r1,v1,r2,v2,r3,v3 = u.reshape([-1,2])
    m1,m2,m3 = masses
    f12, f13, f23 = force(r1-r2), force(r1-r3), force(r2-r3)
    a1,a2,a3 = -m2*f12-m3*f13, m1*f12-m3*f23, m1*f13+m2*f23
    return np.concatenate([v1,a1,v2,a2,v3,a3])

Then the execution I essentially copied your code, just adding some options that are nice for the graph

from scipy import integrate

#Time start, step, and finish point
t0,tf,t_step = 0, 2, 0.1
nsteps = int((tf - t0)/t_step)
u0 = np.concatenate([r1,v1,r2,v2,r3,v3])

solution = integrate.RK45(odesys, t0, u0, tf, first_step=0.2*t_step, max_step=t_step)

# collect data
t_values = [t0]
u_values = [u0]

#The loop for running the Runge-Kutta method over some time period.
for step in range(nsteps):
    solution.step()
    t_values.append(solution.t)
    u_values.append(solution.y)
    # break loop after modeling is finished
    if solution.status == 'finished':
        break

There were no errors reported, and the solution plots as

My scipy version is 1.4.1.

The plot is obtained via

u = np.asarray(u_values).T
# x1,y1 = u[0], u[1],  x2,y2 = u[4],u[5]
plt.plot(u[0],u[1],'-o',lw=1, ms=3, label="body 1")
plt.plot(u[4],u[5],'-x',lw=1, ms=3, label="body 2")
plt.plot(u[8],u[9],'-s',lw=1, ms=3, label="body 3")

Instead of, e.g., u[4] after the transformation of the view of the data in u_values one could also have used u_values[:][4] using the original result data structure.

With changed data so that the center-of-mass is largely constant, and the first body small to the other two, and an increased gravity constant

masses = [0.01, 1, 1]
r1,v1 = [ 0,0], [-2**0.5,-2**0.5]
r2,v2 = [-1,0], [ 0.5**0.5, 0.5**0.5]
r3,v3 = [ 1,0], [-0.5**0.5,-0.5**0.5]
G = 4

The resulting dynamic is