How to solve numerically a nonlinear ODE (a BVP) with singular points in python?

Question

So, I have to solve numerically the following ODE y''+f(y)*(y')^2 = 0: ODE

originally between [y_i,y_f] where y = y(t) and y(0) = y_i. The main problem is that the function f(y) has a (regular) singular point (couldn't post the image). Also, f(y) is obtained from long numerical calculations, so I have no analytical expression for it. The issue is that the singular point lies between y_i and y_f. So far, I could not find any method which helps me to solve this kind of problem. It doesn't matter if it can be solved as a BVP o IVP, but I need to cross the singularity.

What I have tried:

1. I approximated f(y) as -2/(y-y*), where y* is the singularity, and tried to solved the problem as a IVP using Runge-kutta, odeint and Runge-Kutta-Fehlberg method. But I cannot cross it.

2. I tried to cheat using RKF, as y tends to be constant I manually increased it, forcing it to cross the singular point, but then the solution increase to infinity, so not a valid solution.

3. Same as before, but using the numerically defined function f(y), not its approx.

4. The same but as BVP using the shooting method and solve_bvp from Python.

Answer 1

You can immediately integrate after dividing by y' to get ln(y')+F(y)=c or y'*exp(F(y))=C, F'=f, so that in principle you could compute t(y) via simple quadrature, and then get the solution y(t) via inversion of the function value table.

Your example integrates to y'=C*(y-y*)^2, y(t)=y(0)/(1-C*y(0)*t), which could also lead to a singularity in the solution formula due to the denominator. This means that even if a solution exists, you will have to start quite close to it, else the solution process might have to cross a surface of singular Jacobians or other impossibilities, and will fail to continue at this point.