title: Solve system of coupled differential equations using scipy's solve_bvp
body: <p>I want to solve a boundary value problem consisting of 7 coupled 2nd order differential equations. There are 7 functions, <code>y1(x),...y7(x)</code>, and each of them is described by a differential equation of the form</p>
  <pre><code>d^2yi/dx^2 = -(1/x)*dyi/dx - Li(y1,...,y7) for 0 &lt; a &lt;= x &lt;= b,
  </code></pre>
  <p>where <code>Li</code> is a function that gives a linear combination of <code>y1,...,y7</code>. We have boundary conditions for the first-order derivatives <code>dyi/dx</code> at <code>x=a</code> and for the functions <code>yi</code> at <code>x=b</code>:</p>
  <pre><code>dyi/dx(a) = Ai,
  yi(b) = Bi.
  </code></pre>
  <p>So we can rewrite this as a system of 14 coupled 1st order ODEs:</p>
  <pre><code>dyi/dx = zi,
  dzi/dx = -(1/x)*zi - Li(y1,...,y7),
  
  zi(a) = Ai,
  yi(b) = Bi.
  </code></pre>
  <p>I want to solve this system of equations using the Python function <code>scipy.integrate.solve_bvp</code>. However, I have trouble understanding what exactly should be the input arguments for the function as described in the documentation (<a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html" rel="nofollow noreferrer">https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html</a>).</p>
  <p>The first argument that this function requires is a callable <code>fun(x,y)</code>. As I understand it, the input argument <code>y</code> must be an array consisting of values of <code>yi</code> and <code>zi</code>, and gives as output the values of <code>zi</code> and <code>dzi/dx</code>. So my function would look like this (pseudocode):</p>
  <pre><code>def fun(x,y):
      y1, z1, y2, z2, ..., y7, z7 = y
      return [z1, -(1/x)*z1 - L1(y1,...,y7),
              ...,
              z7, -(1/x)*z7 - L7(y1,...,y7)]
  </code></pre>
  <p>Is that correct?</p>
  <p>Then, the second argument for <code>solve_bvp</code> is a callable <code>bc(ya,yb)</code>, which should evaluate the residuals of the boundary conditions. Here I really have trouble understanding how to define such a function. It is also not clear to me what exactly the arrays <code>ya</code> and <code>yb</code> are and what shape they should have?</p>
  <p>The third argument is <code>x</code>, which is the 'initial mesh' with a shape <code>(m,)</code>. Should <code>x</code> consist only of the points <code>a</code> and <code>b</code>, which is where we know the boundary conditions? Or should it be something else?</p>
  <p>Finally the fourth argument is <code>y</code>, which is the 'initial guess for the function values at the mesh nodes', and has shape <code>(n,m)</code>. Its ith column corresponds with <code>x[i]</code>. I guess that the 1st row corresponds with <code>y1</code>, the 2nd with <code>z1</code>, the 3rd with <code>y2</code>, etc. Is that right? Furthermore, which values should be put here? We could put in the known boundary conditions at <code>x=a</code> and <code>x=b</code>, but we don't know how the function looks like at any other points. Furthermore, how does this <code>y</code> relate to the function <code>bc(ya,yb)</code>? Are the input arguments <code>ya,yb</code> somehow derived from this <code>y</code>?</p>
  <p>Any help with understanding the syntax of <code>solve_bvp</code> and its application in this case would be greatly appreciated.</p>