我試圖使用 integrate.solve_ivp 繪制木星圍繞太陽的軌道（以及拉格朗日點中的 2 個星團），但是當我繪制 x 位置和 y 的圖形時，我得到的是螺旋線，而不是穩定的軌道。誰能幫忙？import numpy as npimport scipy.integrateimport matplotlib.pyplot as plt#takes the position of two masses and outputs the vector difference, and the modulusdef vectorise(x1,y1,z1,x2,y2,z2):    v = np.array([x2-x1,y2-y1,z2-z1])    return v, np.linalg.norm(v)def derivatives(t, y):    G =4*np.pi**2    Mj = 0.001    Ms = 1    #Vij gives the vector pointing from i to j (leading to force on j from i)    Vjs = vectorise(y[3],y[4],y[5],y[0],y[1],y[2])    Vsg = vectorise(y[0],y[1],y[2],y[6],y[7],y[8])    Vjg = vectorise(y[3],y[4],y[5],y[6],y[7],y[8])    Vst = vectorise(y[0],y[1],y[2],y[9],y[10],y[11])    Vjt = vectorise(y[3],y[4],y[5],y[9],y[10],y[11])    return [y[12],y[13],y[14],#first differentials of sun position    y[15],y[16],y[17],#first differentials of Jupiter position    y[18],y[19],y[20],#first differentials of Greek position    y[21],y[22],y[23], #first differentials of Trojan position    -G*Mj*1/(Vjs[1]**3) *Vjs[0][0], #second differentail of y[12] (sun x)    -G*Mj*1/(Vjs[1]**3) *Vjs[0][1], #second differentail of y[13] (sun y)    -G*Mj*1/(Vjs[1]**3) *Vjs[0][2], #second differentail of y[14] (sun z)    G*Ms*1/(Vjs[1]**3) *Vjs[0][0], #second differentail of y[15] (Jupiter x)    G*Ms*1/(Vjs[1]**3) *Vjs[0][1], #second differentail of y[16] (Jupiter y)    G*Ms*1/(Vjs[1]**3) *Vjs[0][2], #second differentail of y[17] (Jupiter z)