I am trying to numerically solve equation (6) of [Lakhina 2021][1] in Python. The equation is

$$\frac{1}{2}\left(\frac{d \phi}{d\xi}\right)^2 + S(\phi, M) = 0\, .$$

The Sagdeev potential expression is given by (7). 

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

What I want is to reproduce the potential profiles in Fig. 3 of [Lakhina 2021][1]. 

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

The boundary conditions given in the paper are:

$\phi(0)_{M = 2.55} = 0.023$, $\phi(10)_{M = 2.55} = 0$ <br/>
$\phi(0)_{M = 2.57} = 0.037$, $\phi(12)_{M = 2.57} = 0$ <br/>
$\phi(0)_{M = 2.55} = 0.046$, $\phi(14)_{M = 2.59} = 0$


In the code below, I first define a function for the first-order differential equation. Then, set boundary conditions for each mach number, $M$, and finally, I use `odeint` from the `scipy.integrate` module in Python to solve the boundary value problem. The plot of the solutions is shown in the last figure.

Here is my attempt, the Python code:

```python
## 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 
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_M255 = 0.023 #For M = 2.55
phi0_M257 = 0.037 #For M = 2.57
phi0_M259 = 0.046 #For M = 2.59



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()
```

Output:

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


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

### References

1. Lakhina, G. S., Singh, S. V., & Rubia, R. (2021). [A mechanism for electrostatic solitary waves observed in the reconnection jet region of the Earth’s magnetotail.][1] Advances in Space Research.


  [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