I am solving a ordinary differential equation with scipy‘s solve_bvp to get a family of eigenstates phi_n and eigenvalues omega_n. The equation looks something like
y'' = -(omega_n - f(x)) ** 2 * y,
and I want it to obey the boundary conditions
y(0)=y(1)=0.
f(x) is a given function that depends on some physical parameters.

The problem is that for certain ‘critical’ f(x) the system starts behaving very weirdly, and a bit unstable. Near these instabilities, solve_bvp either doesn’t converge, or the omega_n go to wrong eigenvalues. For example if the eigenvalues are pi, 2pi, 3pi, … and I input as a omega_n pi-0.3, it might converge to 2pi. This is very relevant, as having twice repeated solutions screws over some physical quantities I want to calculate from there.

Also, the solution and eigenvalue guesses are very good. Their either calculated by hand or calculated from previous solve_bvp solutions.

My questions is:

1. What is the integration method solve_bvp uses? Is there a method that might work better here?

I can’t find this in solve_bvp documentation, and I fear I can’t change it.

2. What are things that usually affect the convergence of the parameters? Increasing max_nodes changes the error from ‘repeated eigenvalue’/’solution did not converge’ to MemoryError. Decreasing tol to have better accuracy doesn’t seem to help much either.

I can gladly post the MWE, however I do not find it necessary/useful since:

1. It is about 400-500 lines long,

2. I am asking about the general behaviour of solve_bvp