https://qiita.com/trami/items/a1bdec427fc0420e6f19
Using this site as a reference, I made a simulation of a triple pendulum as shown below, but it can only be calculated up to 0.7 seconds.
In the following equation of motion, the angle is calculated by obtaining the angular velocity, θ, from the first derivative of the angular velocity.
However, no matter how much you investigate, you cannot find the cause. I would appreciate it if you could tell me what caused it.
Thanks for your consideration.

from numpy import sin, cos
from scipy.integrate import solve_ivp
import numpy as np
import csv
from fractions import Fraction as fr
Pi = 3.14 # Pi
G = 9.806 #Gravity acceleration
f = open (result.csv, 'a', newline = "")
writer = csv.writer (f)
th1_DEG = 90
th2_DEG = 90
th3_DEG = 90
om1_DEG = 0
om2_DEG = 0
om3_DEG = 0
L1 = 1
L2 = 1
L3 = 1
M1 = 1
M2 = 1
M3 = 1
th1 = th1_DEG * Pi/180 #convert initial angle to radians
th2 = th2_DEG * Pi/180
th3 = th3_DEG * Pi/180
om1 = om1_DEG * Pi/180
om2 = om2_DEG * Pi/180
om3 = om3_DEG * Pi/180
def func (t, state): # equation of motion
    dydx = np.zeros_like (state)
    dydx [0] = state [1] #dydx [0] = om1
    delta1 = state [0] -state [2] #Delta 1delta2 = state [2] -state [4] #Delta 2
    delta3 = state [4] -state [0] # delta3
    dydx [1] = ((M3 * cos (delta2) -fr (1,4) * (M2 + 4 * M3)) * (state [3] * state [3] * L3 * M3 * sin (delta3)- state [3] * state [3] * (L2 * (M2 + 2 * M3) *
    sin (delta1) + 2 * L2 * M3 * sin (delta2) * cos (delta3)) + state [1] * state [1] * 2 * L1 * M3 * sin (delta3) * cos (delta3) + (2 * M3 * cos (delta3) *
    sin (state [4])-(M1 + 2 * M2 + 2 * M3) * sin (state [0])) * G))/(L1 * ((-M3 * sin (delta2) * sin (delta2) -fr (1,4) * M2) * (2 * M3 * sin (delta3) *
    sin (delta3) + fr (1,2) * M1 + 2 * M2) + fr (1,2) * (M2 * cos (delta1) + 2 * M3 * cos (delta1) -2 * M3 * cos (delta2 ) * cos (delta3)) * (M2 * cos (delta1) +
    2 * M3 * cos (delta1) -2 * M3 * cos (delta2) * cos (delta3))))
    dydx [2] = state [3] #dydx [2] = om2
    dydx [3] = (state [1] * state [1] * L1 * (2 * M3 * M3 * sin (delta3) * cos (delta2) * cos (delta3) * cos (delta3) + M3 * (M2 + 2 * M3) * sin (delta1) *
    cos (delta3) * cos (delta3) + fr (1,2) * M3 * (M2 + 4 * M3) * sin (delta3) * cos (delta3) -fr (1,2) * M3 * (M1 + 4 * M2 + 4 * M3) -fr (1,4) * (M2 + 2 * M3) *
    (M1 + 4 * M2 + 4 * M3) * sin (delta1) -2 * M3 * M3 * sin (delta3) * cos (delta2) * cos (delta2) * cos (delta3)) + state [3] * state [3] * L2 * (fr (1,2) *
    M3 * (M1 + 4 * M2 + 4 * M3) * sin (delta1)))/(L2 * (fr (1,2) * M3 * (M1 + 4 * M2 + 4 * M3) * cos (delta2) * cos (delta2) + fr (1,2) * M3 * (M2 + 4 * M3) *
    cos (delta3) * cos (delta3) + fr (1,4) * (M2 + 2 * M3) * (M2 + 4 * M3) * cos (delta1) + M3 * M3 * (cos (delta2) ** 3 ) * cos (delta3) -2 * M3 * M3 * cos (delta2) *
    cos (delta2) * cos (delta3) * cos (delta3) -fr (1,8) * (M2 + 4 * M3) * (M1 + 4 * M2 + 4 * M3) -M3 * (M2 + 2 * M3 ) * cos (delta1) * cos (delta2) * cos (delta2)-
    fr (1,4) * M3 * (M2 + 4 * M3) * cos (delta2) * cos (delta3)))
    dydx [4] = state [5] #dydx [4] = om3
    dydx [5] = (state [1] * state [1] * ((fr (1,2) * (M2 + 4 * M3) * cos (delta3)-(M2 + 2 * M3) * cos (delta1) * cos (delta2)) * (fr (1,2) * L1 * (M1 + 4 * M2 + 4 *
    M3) * sin (delta3)) + ((M2 + 2 * M3) * cos (delta1) * cos (delta3) -fr (1,2) * (M1 + 4 * M2 + 4 * M3) * cos (delta2 )) * (M2 + 2 * M3) * sin (delta2))-
    state [3] * state [3] * L2 * ((fr (1,2) * (M2 + 4 * M3) * cos (delta3)-(M2 + 2 * M3) * cos (delta1) * cos (delta2 )) * ((M2 + 2 * M3) * sin (delta1) *
    cos (delta3) + fr (1,2) * (M1 + 4 * M2 + 4 * M3) * sin (delta2))-((M2 + 2 * M3) * cos (delta1) * cos (delta3) -fr (1,2) * (M1 + 4 * M2 + 4 * M3) * cos (delta2))
    * (M2 + 2 * M3) * sin (delta2) * cos (delta1)) + state [5] * state [5] * L3 * (((M2 + 2 * M3) * cos (delta1) * cos (delta3 ) -fr (1,2) * (M1 + 4 * M2 + 4 * M3) *
    cos (delta2)) * M3 * sin (delta2) * cos (delta3) + (fr (1,2) * (M2 + 4 * M3) * cos (delta3)-(M2 + 2 * M3) * cos (delta1 ) * cos (delta2)) * M3 * sin (delta3) *
    cos (delta3))-G * ((fr (1,2) * (M2 + 4 * M3) * cos (delta3)-(M2 + 2 * M3) * cos (delta1) * cos (delta2)) * ( (M1 + 2 * M2 + 2 * M3) * sin (state [0]) *
    cos (delta3) -fr (1,2) * (M1 + 4 * M2 + 4 * M3) * sin (state [4]))-((M2 + 2 * M3) * cos (delta1) * cos (delta3 ) -fr (1,2) * (M1 + 4 * M2 + 4 * M3) * cos (delta2)) *
    (M2 + 2 * M3) * cos (state [0]) * sin (delta2)))/(L3 * ((fr (1,2) * (M2 + 4 * M3) * cos (delta3)-(M2 + 2 * M3) * cos (delta1) * cos (delta2)) * (M3 *
    cos (delta3) * cos (delta3) -fr (1,4) * (M1 + 4 * M2 + 4 * M3))-((M2 + 2 * M3) * cos (delta1) * cos (delta3) -fr (1,2) * (M1 + 4 * M2 + 4 * M3) *
    cos (delta2)) * (M3 * cos (delta2) * cos (delta3) -fr (1,2) * (M2 + 2 * M3) * cos (delta1))))
    return dydx
t_span = [0.0,10] #Time generation
dt = 0.01
t = np.arange (t_span [0],

t_span [1],

dt)
writer.writerow ((float (dt), float (th1_DEG), float (th2_DEG), float (th3_DEG), float (om1_DEG), float (om2_DEG), float (om3_DEG), float (L1), float (L2), float (L3), float (M1), float (M2), float (M3)])state = [th1, om1, th2, om2, th3, om3]
print (len (t))
y = solve_ivp (func, t_span, state, t_eval = t) .y
for tt, th1, th2, th3 in zip (t, y [0 ,:],

y [2 ,:],

y [4 ,:]):
    x1 = L1 * sin (th1)
    y1 = -L1 * cos (th1)
    x2 = x1 + L2 * sin (th2)
    y2 = y1-L2 * cos (th2)
    x3 = x2 + L3 * sin (th3)
    y3 = y2-L3 * sin (th3)
    X1 = fr (1,2) * L1 * sin (th1)
    Y1 = -fr (1,2) * L1 * cos (th1)
    X2 = x1 + fr (1,2) * L2 * sin (th2)
    Y2 = y1-fr (1,2) * L2 * cos (th2)
    X3 = x2 + fr (1,2) * L3 * sin (th3)
    Y3 = y2-fr (1,2) * L3 * cos (th3)
    Row_list = [float (tt), float (x1), float (y1), float (X1), float (Y1), float (x2)
    , float (y2), float (X2), float (Y2), float (x3), float (y3), float (X3), float (Y3)]
    writer.writerow (Row_list)
f.close ()

Then

y = solve_ivp (func, t_span, state, t_eval = t) .y

the code here

sol = solve_ivp (func, t_span, state, t_eval = t)
y = sol.y
print (sol.status)

When rewritten like this, -1 was output. When I examine it, it looks like an error called Integration step failed, but I don't know how to solve it.