Issues

• The order of the equation is 2, thus the dimension of the state vector is 2, the value is always y[0], the derivative y[1], there is no y[2], probably a remnant from the translation from Matlab.

• Also in the boundary conditions, there is no ya[2], the derivative value is ya[1], and the function value in the second one is yb[0].

• The initial solution guess has to have the same number 2 of state components.

• Why do you compute two solutions with identical data?

• Remark: It is not necessary to convert the return values into numpy types, the solver checks and converts anyway.

BVP framework with singularity treatment

The ODE is singular at r=0, so one has to treat the first segment in a special way. The mean value theorem gives

(y'(r)-y'(0))/r->y''(0)  for  r->0,

so that in that limit r->0 you get

3*y''(0) = a*y(0) + b*y(0)^3`. 

This allows to define the first arc as

y(r) = y0 + (a*y0 + b*y0^3)*r^2/6
y'(r) = (a*y0 + b*y0^3)*r/3

up to order r³ and r². So if you want precision 1e-9 in y(r), then the first segment should not be longer than 1e-3.

Instead of trying to eliminate y0 from the equations for y(h) and y'(h) to get a condition connecting ya[0] and ya[1], let the solver do also this work and add y0 as parameter to the system. Then the boundary conditions have 3 slots corresponding to the added virtual dimension for the parameter, which can be naturally filled with the equations y(h)=ya[0] and ya[1]=y'(h) and the right boundary condition.

Altogether you can define the system as

h = 1e-3;

def fun(r, y, p):
    return  y[1], -2/r*y[1]+y[0]*(a+b*y[0]*y[0]) 

def bc(ya, yb, p):
    y0, = p
    yh = y0 + y0*(a+b*y0*y0)*h*h/6;
    dyh = y0*(a+b*y0*y0)*h/3
    return ya[0]-yh, ya[1]-dyh, yb[0]-c


x = np.linspace(h, 10, 10)
ya = np.zeros((2, x.size))

sol = solve_bvp(fun, bc, x, ya, p=[1])
print(sol.message,f"y(0)={sol.p[0]}");
plt.plot(sol.x, sol.y[0]);

so that with example parameters a, b, c = -1, 0.2, 3 you get a convergent solver call with y(0)=2.236081087849196 and the resulting plot