我正在尝试用 Python 解决边界值问题。我一直在使用 scipy.integrate.solve_bvp 但它给我的结果是完全错误的。基本上我的代码如下：

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

这个想法是，使用我提供给代码的参数，我应该获得类似于本文中图 [1]（上）的内容：https ://arxiv.org/pdf/0908.2435.pdf ，但是如果我运行我得到的代码完全错误。我可以使用另一个求解器吗？我的代码有问题吗？

谢谢你的时间。问候，

路易斯·P。