Using Velocity Verlet to simulate 3-body motion in 3D

Question

As the title states, I want to use the velocity Verlet algorithm to simulate the trajectories of 3 bodies in 3D. The code is based from this. My problem is the trajectories outputted seem to be small oscillations propagating as straight lines, whereas I should be getting something chaotic like random squiggles. I'm not sure if the problem lies in my implementation of the velocity Verlet algorithm or in the mathematics defining the gravitational force between 3 bodies or whatever else.

import numpy as np 
import matplotlib.pyplot as plt 
import numpy.linalg as lag
from mpl_toolkits.mplot3d import Axes3D
plt.style.use('dark_background')

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

# starting coordinates for suns

r1_initial = np.array([-10e10, 10e10, -11e10])
v1_initial = np.array([3e10, 0, 0])

r2_initial = np.array([0, 0, 0])
v2_initial = np.array([0, 3e10, 0])

r3_initial = np.array([10e10, 10e10, 12e10])
v3_initial = np.array([-3e10, 0, 0])

def force_grav_3bodies(r1,r2,r3,m1,m2,m3,G):
    force_A = -G*m1*m2*(r1-r2)/lag.norm(r1-r2)**3 - G*m1*m3*(r1-r3)/lag.norm(r1-r3)**3
    
    force_B = -G*m2*m1*(r2-r3)/lag.norm(r2-r3)**3 - G*m2*m3*(r2-r1)/lag.norm(r2-r1)**3
    
    force_C = -G*m3*m1*(r3-r1)/lag.norm(r3-r1)**3 - G*m3*m2*(r3-r2)/lag.norm(r3-r2)**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

dt = 0.1
Tmax = 2000
steps = int(Tmax/dt)
G = 6.67e-11

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)

# velocity verlet algorithm 

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

    v1[k] = v1[k-1] + 0.5*(a1[k-1] + a1[k])*dt
    v2[k] = v2[k-1] + 0.5*(a2[k-1] + a2[k])*dt
    v3[k] = v3[k-1] + 0.5*(a3[k-1] + a3[k])*dt
    
    
fig = plt.figure(figsize=(8, 8))
ax = fig.gca(projection='3d')
plt.gca().patch.set_facecolor('black')

plt.plot([i[0] for i in r1], [j[1] for j in r1], [k[2] for k in r1] , '^', color='r', lw = 0.5, markersize = 0.1, alpha=0.5)
plt.plot([i[0] for i in r2], [j[1] for j in r2], [k[2] for k in r2] , '^', color='y', lw = 0.5, markersize = 0.1, alpha=0.5)
plt.plot([i[0] for i in r3], [j[1] for j in r3], [k[2] for k in r3] , '^', color='c', lw = 0.5, markersize = 0.1, alpha=0.5)

ax.scatter(r1[steps-2][0],r1[steps-2][1],r1[steps-2][2],color="r",marker="^",s=100)
ax.scatter(r2[steps-2][0],r2[steps-2][1],r2[steps-2][2],color="y",marker="^",s=100)
ax.scatter(r3[steps-2][0],r3[steps-2][1],r3[steps-2][2],color="c",marker="^",s=100)

ax.scatter(r1[0][0],r1[0][1],r1[0][2],color="r",marker="o",s=100,label="A")
ax.scatter(r2[0][0],r2[0][1],r2[0][2],color="y",marker="o",s=100,label="B")
ax.scatter(r3[0][0],r3[0][1],r3[0][2],color="c",marker="o",s=100,label="C")
plt.legend()
plt.axis('on')

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()