I'm using scipy.integrate.solve_ivp to solve a system of ODE and it works fine. Now I'm introducing the spatial dependency and I'm having some issues.
This is a simplification of my system (diffusion plus reactions):
The time derivative is what fun returns (solve_ivp fun parameter). For the spatial derivative, I use a finite-different formula (by hand). For each position you have this set of equations, so fun returns the time derivatives of each yi at each x.
The code works fine and the result has sense, but I'm struggling with stability. The system is stiff and I use the Radau method. The problem is related to rtol, atol and the space discretisation, dx. If atol (dominates over rtol) is too big, the system diverges, but if it is too small, it also diverges, I have to also reduce dx to make it work. I think that the problem comes from solve_ivp, which is looking for very small dt to accomplish with atol but since dx cannot be affected by dt, somehow it diverges.
Do you know any general stability relation between dx and dt? Or may I use a different approach?