I am striving to solve a coupled nonlinear ordinary differential equation, or actually a second-order differential equation, with the following form

∂π(t)/∂t=−3coth(3t)*π(t)−φ(t)+φ(t)^3

∂φ(t)/∂t=π(t)

I used the function odeint of python's scipy package. Here are the following code

import numpy as np
from scipy.integrate import odeint
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

#function that return dz/dt
def model(z, t):
    phi=z[0]
    pi=z[1]
    dVdphi = phi**3
    h = 3/(np.tanh(3*t))
    dphidt = pi
    dpidt = -phi+dVdphi-h*pi
    dzdt = [dphidt, dpidt]
    return dzdt

phi = np.empty(1000)
pi = np.empty(1000)
for i in range(1,1000) :
    rand = np.random.rand(2)
    #initial condition
    z0 = np.empty(2)
    z0[0]=1.5*rand[0]
    z0[1]=0.4*rand[1]-0.2
    #number of time points
    n=1000
    #time points
    t = np.linspace(1,0.0001,n)
    z = odeint(model, z0, t, mxstep=500000)
    phi[i] = z[n-1][0]
    pi[i] = z[n-1][1]


plt.scatter(phi, pi, label ='positive')
plt.ylabel('pi')
plt.xlabel('phi')
plt.savefig('out.png', dpi=150)
plt.show()

However, whenever I run the code, it keeps showing me Isoda error like

"lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step
       (h = step size)."

Is there anyone who can help me solve this problem? Should I use another tool to solve this differential equation?