Error: _vhstack_dispatcher() takes 1 positional argument but 9 were given

[alexfrol86]

Hi!

I have a system of 10 nonlinear ODE with 10 variables (y[0], y[1], ... y[9]), there are 10 initial conditions for y[i] for x=0 and 4 conditions for x = 6.804.

I try to solve a boundary value problem using scipy.integrate.solve_bvp according to the following algorithm: https://docs.scipy.org/doc/scipy/referen...e_bvp.html
I get the error "_vhstack_dispatcher() takes 1 positional argument but 9 were given" though I tried to use np.dstack instead of np.hstack. Please help to fix the error.

My code is below:

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
from scipy.constants import g

def fun(x, y):
    return np.vstack((ptb / y[4] * y[5] / np.sqrt(y[5]**2 + y[6]**2) * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes) - 1 / y[2]**2 + y[1]**2 / y[2],
                      ptb / y[4] * y[6] / np.sqrt(y[5]**2 + y[6]**2) * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes) - y[0] * y[1] / y[2]),
                     y[0], y[1] / y[2], 
                      -mqb * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes), y[6] * y[1] / y[2] - y[7], -y[5] * 2.0 * y[1] / y[2] + y[6] * y[0] / y[2] - y[8] * 1.0 / y[2],
                      y[5] * (-2.0 / y[2]**3 + y[1]**2 / y[2]**2) + y[6] * (-y[0] * y[1]) / y[2] + y[8] * y[1] / y[2]**2, 0,
                       ptb / y[4]**2 * (np.sqrt(y[5]**2 + y[6]**2)) * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes))
    

#Implement evaluation of the boundary condition residuals:

def bc(ya, yb):
    return np.array([ya[-0.03899746484, 1.05800280714, -0.30648856279, -0.605427, 1, -0.065851, 1.361805, 1.419320, -8.219418e-5, 0.6], yb[3.230021179060866e-17, 0.5519654178809713, 2.221001221001221, np.pi]])
#Define the initial mesh with 20 nodes:

x = np.linspace(0.0, 6.804, 20)

#To obtain solutions, different initial guesses for y are used. They are denoted by subscripts a and b.
y_a = np.zeros((10, x.size))
y_b = np.zeros((10, x.size))
y_b[0] = 3.230021179060866e-17

res_a = solve_bvp(fun, bc, x, y_a)
res_b = solve_bvp(fun, bc, x, y_b)

x_plot = np.linspace(0, 6.804, 100)
y_plot_a = res_a.sol(x_plot)[0]
y_plot_b = res_b.sol(x_plot)[0]
#plt.plot(x_plot, y_plot_a, label='y_a')
#plt.plot(x_plot, y_plot_b, label='y_b')
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.show()