I’m trying to solve a differential equation (actually it’s part of an optimization process, so it’s solved each optimization iteration). The problem is: my u is a parameter that has a physical meaning of extent, so it can only be between 0 and 1. However, sometimes the solver tries to use u < 0.0 or y > 1.0, which doesn’t make (physical) sense and causes the resulting value to be complex.
I followed the advice of saturating the value with u = max(min(u, 1.0), 0.0), but when I do this, the next iteration comes with u = Inf.
I wrote a minimum (but complete) example showing this behavior. The hardcoded values come from the optimizer.
The code is here.
What am I missing? How can I tell OrdinaryDiffEq.jl the domain of my function?
Thank you.