title: Solving a boundary value problem (diffusion-reaction equation) with scipy solve_bvp
body: <p>I'm struggling to solve the following 2nd order boundary value problem:</p>
  
  <pre><code>y'' + 2/x*y' + k**2.0*F(y) = 0
  
  y(x=1)=1,  y'(x=0)=0
  
  F(y) = -y or F(y) = -y*exp(AB*(1-y)/(1+B(1-y))
  </code></pre>
  
  <p>I somehow fail to set the boundary conditions right. I defined the function for <code>F(y)=y</code> and boundary conditions the following way:</p>
  
  <pre><code>def fun(x, y, p):
   k = p[0]
   return np.vstack((y[1], -2.0/x*y[1] + k**2.0*y[0]))
  
  def bc(ya, yb, p):
       return np.array([ya[0], yb[0],ya[1]])
  
  y[0,:] = 1
  y[1,0] = 0
  sol = solve_bvp(fun, bc, x, y, p=[40])
  </code></pre>
  
  <p>The results should I get are definitely wrong, changing the initial conditions doesn't make things better. I think my problem is some how related to the zero gradient boundary condition at x=0. Does anybody know what I'm doing wrong?</p>
  
  <p>EDIT:
  Here a MWE, which should give a constant value of 1 for k=0.01. but for a k=5, the value at x=0 should be approx. 0.06:</p>
  
  <pre><code>def fun(x, y, p):
       k = p[0]
       return np.vstack((y[1], -2.0/x*y[1] + k**2.0*y[0]))
  
  def bc(ya, yb, p):
       return np.array([ya[0], yb[0]-1.0,yb[1]])
  
  x = np.linspace(1e-3, 1, 100)
  y = np.zeros((2, x.size))
  y[0,:] = 1
  
  from scipy.integrate import solve_bvp
  sol = solve_bvp(fun, bc, x, y, p=[1000])
  x_plot = np.linspace(0, 1, 100)
  y_plot = sol.sol(x_plot)[0]
  plt.figure()
  plt.plot(x_plot, y_plot)
  </code></pre>