I am trying to solve equation (6) of Lakhina 2021, numerically, in Python, so that I can reproduce the potential profiles in Fig. 3 of Lakhina 2021. The Sagdeev potential expression is given by (7).

In the code below, I first define a function for the ode. Then, set an initial boundary condition, and finally, I use odeint from the scipy.integrate module in Python to solve the ode. The plot of the solution is shown in the last figure.

Here is my attempt:

##Importing standard modules
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt 


##Reconnection jet plasma parameters
n1 = 0.74 
n2 = 0.26 
sig1 = 0.11 
sig2 = 0.07 
U1 = -1.72
U2 = 1.82 


#Function for Sagdeev potential equation (fast)
def S(phi, M):
    s = (1 - np.exp(phi)) + n1/(6*np.sqrt(3*sig1))*((M - U1 + np.sqrt(3*sig1))**3 -
                                                  ((M - U1 + np.sqrt(3*sig1))**2 - 2*phi)**1.5 -
                                                  (M -  U1 - np.sqrt(3*sig1))**3 + 
                                                  ((M - U1 - np.sqrt(3*sig1))**2 - 2*phi)**1.5) + n2/(6*np.sqrt(3*sig2))*(
                                                  (M - U2 + np.sqrt(3*sig2))**3 -
                                                  ((M - U2 + np.sqrt(3*sig2))**2 - 2*phi)**1.5 -
                                                  (M - U2 - np.sqrt(3*sig2))**3 +
                                                  ((M - U2 - np.sqrt(3*sig2))**2 - 2*phi)**1.5) 
    return s

##Solving the ode

def model(phi, zeta, M):

S = (1 - np.exp(phi)) + n1/(6*np.sqrt(3*sig1))*((M - U1 + np.sqrt(3*sig1))**3 -
                                              ((M - U1 + np.sqrt(3*sig1))**2 - 2*phi)**1.5 -
                                              (M -  U1 - np.sqrt(3*sig1))**3 + 
                                              ((M - U1 - np.sqrt(3*sig1))**2 - 2*phi)**1.5) + n2/(6*np.sqrt(3*sig2))*(
                                              (M - U2 + np.sqrt(3*sig2))**3 -
                                              ((M - U2 + np.sqrt(3*sig2))**2 - 2*phi)**1.5 -
                                              (M - U2 - np.sqrt(3*sig2))**3 +
                                              ((M - U2 - np.sqrt(3*sig2))**2 - 2*phi)**1.5)  
dphi_dzeta = -np.sqrt(-2*S)


return dphi_dzeta


#Boundary conditions
phi0 = 0.023



phi_array = np.linspace(-0.01, 0.06, 1000)
zeta_array = np.linspace(-16, 16, 1000)

Phi = odeint(model, phi0, zeta_array, args = (2.57,))

##Plotting

plt.figure(2)
plt.axhline(0, color = 'k', lw = 1)
plt.axvline(0, color = 'k', lw = 1)
plt.plot(zeta_array, Phi, label = "M = 2.55")
plt.xlabel("$\zeta$")
plt.ylabel("S($\phi$, M)")
plt.legend()

Ouput:

May you please assist? I am really not sure where I am going wrong. Thank you in advance.