Which solvers for BVP in python are the best? Is there something better that scipy.integrate.solve_bvp?

I am trying to solve a boundary value problem with Python. I have been using scipy.integrate.solve_bvp but the result that it is giving me is completely wrong. Basically my code is as follows:

from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import solve_bvp

r_in = 0.0
r_end = 500.
l = 0.
rho_c = 1.
m = 1.

def system_DE(r,y,p_self,l=l,rho_c=rho_c):
R = y[0]
ph = y[1]
a = y[2]
al = y[3]

om = p_self[0]
u = p_self[1]

dR_dr = ph
da_dr = a*(((2.*l+1.)*r/2.)*(-om**2.*a**2.*R**2./al**2.+ph**2.+l*(l+1.)*a**2.*R**2./r**2.+m**2.*a**2.*R**2.)-(a**2.-1.)/(2.*r))
dal_dr = al*(da_dr/a-l*(l+1.)*(2.*l+1.)*a**2.*R**2./r-(2.*l+1.)*m**2.*a**2.*r*R**2.+(a**2.-1.)/r)
dph_dr = -2.*ph/r+(l*(l+1.)*(2.*l+1.)*a**2.*R**2./r+(2.*l+1.)*m**2*a**2*r*R**2-(a**2-1)/r)*ph-om**2.*a**2.*R/al**2.+l*(l+1.)*a**2.*R/r**2.+m**2.*a**2.*R

return [dR_dr,dph_dr,da_dr,dal_dr]

def bc(ya,yb,p_sch,l=l,rho_c=rho_c,r_in=r_in,r_end = r_end,m=m):

om = p_self[0]
u = p_self[1]

if l==0.:
return np.array([ya[0]-rho_c,ya[1]-rho_c*r_in,yb[0]*((m**2.-om**2.)**(1./2.)+1./(r_end))+yb[1],ya[2]-1.,ya[3]-u,yb[3]*yb[2]-1.])
else:
return np.array([rho_c*r**l,rho_c*l*r**(l-1.),1.,u])

n_r = 1000
r = np.linspace(r_in,r_end,n_r)
y = np.zeros((4,r.size))
y[0,0] = rho_c
y[0,1] = rho_c

y[2,:] = 1.

p_self=[0.8,0.85] #[om,u]

sol = solve_bvp(system_DE,bc, r, y, p=p_self, tol=0.01,max_nodes=100000)

x_plot = np.linspace(r_in, 200  , 5000)
L = sol.sol(x_plot)[0]
Ph = sol.sol(x_plot)[1]
a = sol.sol(x_plot)[2]
alpha = sol.sol(x_plot)[3]

plt.subplot(211)
plt.plot(x_plot, L)
plt.xlabel("x")
plt.ylabel("y")
plt.subplot(212)
plt.plot(x_plot, a)
plt.plot(x_plot, alpha)
plt.xlabel("x")
plt.ylabel("y")
plt.show()


The idea is that with the parameters that I am giving to the code I should obtain something similar to the figure [1] (up) on this article: https://arxiv.org/pdf/0908.2435.pdf, however if I run the code I am obtaining something completely wrong. Is there another solver that I could use? Is there something wrong with my code?

• Check your code, I am using Spyder and I am seeing some nan values for Ph, a and alpha. Are you sure your implementation is correct? – David Oct 14 '18 at 2:34
• Why do you use the global p_self in bc when the passed parameter vector is p_sch? You should consult the fields "sol.success" and "sol.message", there might be something about a "singular Jacobian". – Lutz Lehmann Nov 9 '18 at 9:45