How to solve the integral-like energy equation with Sagdeev potential numerically in Python?

Quite hastily, I admit, I wrote the code below and generated the solutions. I did not bother to center them.

One issue you have with your code is that you are trying to integrate

$$\frac{d\phi}{d\xi} \, =\, \sqrt{- \, 2 \, S(\phi, M)}$$

and the square root may get in the way, because you have to choose whether it is the one above or this one, especially when the derivative becomes zero (and that's why you are not seeing the maximum when you integrate the problem):

$$\frac{d\phi}{d\xi} \, =\, - \,\sqrt{-\,2 \, S(\phi, M)}$$

Because I do not want to deal with this issue, I integrated the Hamiltonian equations instead

\begin{align} &\frac{d\phi}{d\xi} = p\\ &\frac{dp}{d\xi} = - \,\frac{\partial S}{\partial \phi}(\phi, M)\\ \end{align}

In addition to that, I used a python solver for initial conditions (I think that was another issue in your code), so I had to integrate forward in time $[0, \, 16]$ and backward in time $[-16, \, 0]$ and plot both solutions to get the solution from $[-16, \, 16]$

Finally, since the system is time-invariant, if you have a solution $\phi(\xi)$ then $\phi(\xi-\xi_0)$ is also a solution. So I assume that in the paper, the solutions were shifted to have their peaks at $\xi = 0$.

import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def _grad_U(phi, A, a): return - 3 * np.sqrt( A**2 - 2*phi ).dot(a) + np.exp(phi) def f(t, y, A, a): # a = [a1, -a1, a2, -a2] return np.array([ y[1], _grad_U(y[0], A, a) ]) def initialize(phi, A, a): S = (A**3 - np.sqrt( A**2 - 2*phi )**3).dot(a) + 1 - np.exp(phi) return np.sqrt(-2*S) def set_parameters(M, U1, U2, sigma1, sigma2, n1, n2): sqrt_3sigma1 = np.sqrt(3*sigma1) sqrt_3sigma2 = np.sqrt(3*sigma2) A1 = M - U1 + sqrt_3sigma1 A1_ = M - U1 - sqrt_3sigma1 A2 = M - U2 + sqrt_3sigma2 A2_ = M - U2 - sqrt_3sigma2 a1 = n1 / (6*sqrt_3sigma1) a2 = n2 / (6*sqrt_3sigma2) return np.array([A1, A1_, A2, A2_]), np.array([a1, -a1, a2, -a2]) def solve_equation(phi_in, M, U1, U2, sigma1, sigma2, n1, n2, zeta_left, zeta_right, n_zetas): A, a = set_parameters(M, U1, U2, sigma1, sigma2, n1, n2) v_in = initialize(phi_in, A, a) zeta_array = np.linspace(0, zeta_right, n_zetas) y_in = np.array([phi_in, v_in]) sol = solve_ivp(fun = lambda t, y : f(t, y, A, a), t_span=[0, zeta_right], y0=y_in, t_eval=zeta_array, method='Radau') zeta_array = np.linspace(0, zeta_left, n_zetas) sol_ = solve_ivp(fun = lambda t, y : f(t, y, A, a), t_span=[0, zeta_left], y0=y_in, t_eval=zeta_array, method='Radau') return sol, sol_ ## Reconnection jet plasma parameters n1 = 0.74 n2 = 0.26 sig1 = 0.11 sig2 = 0.07 U1 = -1.72 U2 = 1.82 # 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 zeta_left = -16 zeta_right = 16 n_zetas = 5000 sol_255, sol_255_ = solve_equation(phi0_M255, 2.55, U1, U2, sig1, sig2, n1, n2, zeta_left, zeta_right, n_zetas) sol_257, sol_257_ = solve_equation(phi0_M257, 2.57, U1, U2, sig1, sig2, n1, n2, zeta_left, zeta_right, n_zetas) sol_259, sol_259_ = solve_equation(phi0_M259, 2.59, U1, U2, sig1, sig2, n1, n2, zeta_left, zeta_right, n_zetas) #phi_array = np.linspace(-0.01, 0.06, 1000) #zeta_array = np.linspace(-16, 16, 1000) ## Plotting plt.figure(2) plt.axhline(0, color = 'k', lw = 1) plt.axvline(0, color = 'k', lw = 1) plt.plot(sol_255.t, sol_255.y[0,:], color = 'b', label = "M = 2.55") plt.plot(sol_255_.t, sol_255_.y[0,:], color = 'b') plt.plot(sol_257.t, sol_257.y[0,:], color = 'g', label = "M = 2.57") plt.plot(sol_257_.t, sol_257_.y[0,:], color = 'g') plt.plot(sol_259.t, sol_259.y[0,:], color = 'r', label = "M = 2.59") plt.plot(sol_259_.t, sol_259_.y[0,:], color = 'r') plt.xlabel("$\zeta$") plt.ylabel("$\phi$") plt.legend() 

Edit. The first function, _grad_U() calculates $ - \, \frac{\partial S}{\partial \phi}(\phi, M)$, taking as input $\phi$, as well as the parameters featured in the formula of $S$. I wrote it in a vectorized form, that is why it looks shorter (it is easier to keep track of things and make corrections if necessary, plus numpy is faster with vectorized quantities) The second one, f(), takes in $y = (\phi, \, p)$ and outputs $\Big(p, \, - \, \frac{\partial S}{\partial \phi}\partial(\phi, M)\Big )$. Additional parameters are the time t, which is there only because solve_ivp() needs it, and the rest are the parameters for the calculation of the potential. Basically, f() is the right hand side of the differential equations I wrote above. The third function, initialize() takes in the initial state $\phi(0)$ and calculates the corresponding initial momentum $p(0)$ so that they satisfy the total energy equation: $$\frac{1}{2}p(0)^2 \, + \, S(\phi(0), M) \, =\, 0 $$ That is why, the latter function calculates $$p(0) = \sqrt{2 \, S(\phi(0), M)}$$
The forth function, set_parameters() takes in the original parameters that are included in the calculation of the potential $S$ and reorganzies them into the new sets of more convenient parameters A and a, that are fed to the two functions before it.