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

[![enter image description here][2]][2][![enter image description here][3]][3]


[![enter image description here][4]][4]


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:

[![enter image description here][5]][5]


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


  [1]: https://www.sciencedirect.com/science/article/pii/S0273117721003112
  [2]: https://i.sstatic.net/AKLQL.png
  [3]: https://i.sstatic.net/j8Zrk.png
  [4]: https://i.sstatic.net/VPMR6.png
  [5]: https://i.sstatic.net/l499T.png