Rather than implementing by hand a BVP method, one can choose to use SciPy's own BVP solver. The solver is named solve_bvp and takes in input, in its minimal version, four arguments:
A vector function, called ode in the code below, that returns the righ hand side of the first-order ODE system (if the problem has higher order, the order must e brought to 1 with the usual trick);
A vector function, called bc in the code below, with all conditions that must be satisfied. Usually these are just the bounday conditions;
A vector, called t in the code below, for the independent variable (often this is the time variable);
A matrix, called x in the code below, which represents the starting point of the algorithm for the evolutions of the dependent variables.
