Scipyのodeintでエラーが出てしまいます。

前提・実現したいこと

Scipyのodeintを使って二次元ハミルトン系を解析したい

発生している問題・エラーメッセージ

RuntimeWarning: invalid value encountered in double_scalars
ODEintWarning: Excess accuracy requested (tolerances too small). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

該当のソースコード

python
1m = 8.0
2D = 10.0
3z = 2.3
4l = 1.95
5E = 0.10
6
7def model(f,t):
8
9    p1 =f[0]
10    q1 =f[1]
11    q2 =f[2]
12
13    dp1dt = -8*(2*(q1**3)-q1)*np.exp(-z*l*q2)
14    dq1dt = (1/m)*p1
15    dq2dt = ((2/m)*(E-(1/2*m)*(p1**2)-4*((q1**4)-(q1**2))*np.exp(-z*l*q2)+D*(1-np.exp(-l*q2))**2))**(1/2)
16    dfdt = [dp1dt,dq1dt,dq2dt]
17
18    return dfdt
19
20#time
21n1 = 101
22t = np.linspace(0,2,n1)
23p1 = np.empty_like(t)
24q1 = np.empty_like(t)
25q2 = np.empty_like(t)
26
27#initial condition
28n2 = 10
29ini = np.linspace(-1.0,1.0,n2)
30
31for h in range(0,n2):
32    p1_0 = ini[h]
33
34    for i in range(0,n2):
35        q1_0 = ini[i]
36
37        for j in range(0,n2):
38            q2_0 = ini[j]
39
40            f0 =[p1_0,q1_0,q2_0]
41            p1[0] = f0[0]
42            q1[0] = f0[1]
43            q2[0] = f0[2]
44
45            for k in range(1,n1):
46                tspan = [t[k-1],t[k]]
47                f = odeint(model,f0,tspan)
48                p1[k] = f[1][0]
49                q1[k] = f[1][1]
50                q2[k] = f[1][2]
51                f0 = f[1]