3D 3-body motion simulation

Question

I am trying to use the velocity Verlet to simulate the 3-body problem. I am 95% sure that I gravitational force formula and the velocity Verlet algorithm is correct (that is I didn't seem to have any problems when testing them), thus I'm not sure why my plot is incorrect as I'm expecting chaotic lines. Could it be due to the logic in the design/how I implemented the functions for gravitational force and verlet? I'm not very experienced with Python (programming in general) so I really can't grasp why the code won't produce my desired plots.

import numpy as np 
import matplotlib.pyplot as plt 
import numpy.linalg as lag
from mpl_toolkits.mplot3d import Axes3D


# masses of stars
m1 = 10e30
m2 = 20e30
m3 = 30e30

# starting coordinates for stars

r1_initial = np.array([-10e10, 10e10, -10e10])
v1_initial = np.array([5e9, 0, 0])

r2_initial = np.array([0, 0, 0])
v2_initial = np.array([0, 5e9, 0])

r3_initial = np.array([10e10, 10e10, 10e10])
v3_initial = np.array([-5e9, 0, 0])

G = 6.67e-11

def force_grav_3bodies(r1,r2,r3,m1,m2,m3,G):
    force_A = -G*m2*m1*(r2-r1)/lag.norm(r2-r1)**3 - G*m3*m1*(r3-r1)/lag.norm(r3-r1)**3
    
    force_B = -G*m3*m2*(r3-r2)/lag.norm(r3-r2)**3 - G*m1*m2*(r1-r2)/lag.norm(r1-r2)**3
    
    force_C = -G*m1*m3*(r1-r3)/lag.norm(r1-r3)**3 - G*m2*m3*(r2-r3)/lag.norm(r2-r3)**3
    return force_A, force_B, force_C
    
def accel_3sun(r1,r2,r3,m1,m2,m3,G):
    fg3bA, fg3bB, fg3bC = force_grav_3bodies(r1,r2,r3,m1,m2,m3,G) 
    accel_sunA, accel_sunB, accel_sunC = fg3bA/m1, fg3bB/m2, fg3bC/m3
    return accel_sunA, accel_sunB, accel_sunC

# velocity verlet algorithm 

def vverlet_3D_3bodies(Tmax, dt, G):
    steps = int(Tmax/dt)
    
    r1 = np.zeros([steps,3])
    v1 = np.zeros([steps,3])
    a1 = np.zeros([steps,3])

    r2 = np.zeros([steps,3])
    v2 = np.zeros([steps,3])
    a2 = np.zeros([steps,3])

    r3 = np.zeros([steps,3])
    v3 = np.zeros([steps,3])
    a3 = np.zeros([steps,3])
    
    # initial conditions
    
    r1[0], r2[0], r3[0] = r1_initial, r2_initial, r3_initial

    v1[0], v2[0], v3[0] = v1_initial, v2_initial, v3_initial
    
    a1[0], a2[0], a3[0] = accel_3sun(r1[0],r2[0],r3[0],m1,m2,m3,G)
    
    for k in range(steps-1):
    
        r1[k+1] = r1[k] + v1[k]*dt + 0.5*a1[k]*dt*dt
        r2[k+1] = r2[k] + v2[k]*dt + 0.5*a2[k]*dt*dt
        r3[k+1] = r3[k] + v3[k]*dt + 0.5*a3[k]*dt*dt
        
        a1[k+1], a2[k+1], a3[k+1] = accel_3sun(r1[k+1],r2[k+1],r3[k+1],m1,m2,m3,G)

        v1[k+1] = v1[k] + 0.5*(a1[k] + a1[k+1])*dt
        v2[k+1] = v2[k] + 0.5*(a2[k] + a2[k+1])*dt
        v3[k+1] = v3[k] + 0.5*(a3[k] + a3[k+1])*dt
    
    return r1, r2, r3

r1_sol, r2_sol, r3_sol = vverlet_3D_3bodies(8e2,1e-1,G)
    
fig = plt.figure(figsize=(16,16))
ax=plt.axes(projection="3d")
plt.gca().patch.set_facecolor('black')

plt.plot(r1_sol[0:,0], r1_sol[0:,1], r1_sol[0:,2], color='r')
plt.plot(r2_sol[0:,0], r2_sol[0:,1], r2_sol[0:,2], color='y')
plt.plot(r3_sol[0:,0], r3_sol[0:,1], r3_sol[0:,2], color='c')

ax.scatter(r1_sol[0][0],r1_sol[0][1],r1_sol[0][2],color="r",marker="o",s=100,label="A")
ax.scatter(r2_sol[0][0],r2_sol[0][1],r2_sol[0][2],color="y",marker="o",s=100,label="B")
ax.scatter(r3_sol[0][0],r3_sol[0][1],r3_sol[0][2],color="c",marker="o",s=100,label="C")
plt.legend()

ax.w_xaxis.set_pane_color((0.0, 0.0, 0.0, 1.0))
ax.w_yaxis.set_pane_color((0.0, 0.0, 0.0, 1.0))
ax.w_zaxis.set_pane_color((0.0, 0.0, 0.0, 1.0))

plt.show()