我想求解一个由7个耦合的二阶微分方程组成的边值问题，其中有7个函数y1(x),...y7(x)，每个函数都由以下形式的微分方程来描述：

d^2yi/dx^2 = -(1/x)*dyi/dx - Li(y1,...,y7) for 0 < a <= x <= b,

其中Li是给出y1,...,y7的线性组合的函数。我们有dyi/dx在x=a处的一阶导数和yi在x=b处的函数的边界条件：

dyi/dx(a) = Ai,
yi(b) = Bi.

因此，我们可以将其改写为14个耦合的一阶常微分方程的系统：

dyi/dx = zi,
dzi/dx = -(1/x)*zi - Li(y1,...,y7),

zi(a) = Ai,
yi(b) = Bi.

我想使用Python函数scipy.integrate.solve_bvp来求解这个方程组，但是我很难理解文档（https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html）中描述的函数的输入参数到底应该是什么。
这个函数需要的第一个参数是一个可调用的fun(x,y)。据我所知，输入参数y必须是一个由yi和zi的值组成的数组，并给出zi和dzi/dx的值作为输出。因此，我的函数看起来像这样（伪代码）：

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

对不对？
然后，solve_bvp的第二个参数是一个可调用的bc(ya,yb)，它应该计算边界条件的残差。在这里，我真的很难理解如何定义这样一个函数。我也不清楚数组ya和yb到底是什么，它们应该有什么形状？
第三个参数是x，它是形状为(m,)的“初始网格”。x应该只由点a和b组成，这是我们知道边界条件的地方？还是应该是其他的东西？
最后第四个参数是y，它是“网格节点处函数值的初始猜测”，形状为(n,m)。它的第i列对应于x[i]。我猜第一行对应于y1，第二行对应于z1，第三行对应于y2，以此类推。对吗？此外，我们可以在x=a和x=b处输入已知的边界条件，但是我们不知道函数在其他点上是什么样子的。另外，这个y和函数bc(ya,yb)有什么关系？输入参数ya,yb是不是以某种方式从这个y导出的？
对于理解solve_bvp的语法及其在这种情况下的应用的任何帮助都将非常感谢。