I’m trying to solve a problem of filling and emptying of a retention behind a dam. What I’ve got are: 1) a function that relates the elevation of the water in the retention with the volume of the water, V(h); 2) a function that relates the discharge through the dam outlet with the water elevation, q(h); and 3) a function of water inflow given the duration of rain, Vr(tr).
For the case when I ignore the third function, i.e. I assume that V0 appears instantly in the reservoir, I managed to get the solution as:

def dtdh(h, y):
    t = y[0]
    dVdh = misc.derivative(V_h_fun, h, dx=0.1, args=(aCV, bCV, cCV))
    return (- dVdh * 1e3 / Q_h_fun(h, aCQ, bCQ, cCQ)) # V_h_fun relates h[m] to V[dam^3], so 1e3 is to get [m^3]
                                                      # as Q_h_fun gives q[m^3/s] from h[m]

elevs = np.linspace(235,231,100) # initial elev. 235, to final elev. 231 m
sol = integrate.solve_ivp(dtdh, t_span=(235,231), y0=[0], t_eval=elevs) # y0 = t0

V_h_fun and Q_h_fun are the first two functions I mentioned previously, and aCV, bCV, cCV, aCQ, bCQ, cCQ are the coefficients from fitting the curves to some data.
This works great, it gives me the time required to empty the reservoir, given an initial elevation, which is indirectly given through the elevs and t_span.

When I try to expand this to also consider the increase of the volume in the retention with time (due to rain), then the new diff. eq. looks like:

def dtdh2(h, y):
    t = y[0]
    dVdh = misc.derivative(V_h_fun, h, dx=0.1, args=(aCV, bCV, cCV))
    dVrdt = misc.derivative(initial_volume, t, dx=0.1, args=(V0, t_rain))
    dVrdt.clip(min=0, out=dVrdt) # the derivative goes to a large negative value before going back to 0
    return (- dVdh * 1e3 / (Q_h_fun(h, aCQ, bCQ, cCQ)  - dVrdt))

So, I only subtracted the inflow due to the rain, from the outflow. The inflow function is defined as:

def initial_volume(t, V0, t_rain):
    # V0 in [m3]
    # t_rain in [s]
    if type(t) != np.ndarray:
        t = np.array([t])
    
    f = np.zeros_like(t)
    c1 = (t < 0)
    c2 = (t >= 0) & (t <= t_rain)
    c3 = (t > t_rain)
    f[c1] = 0
    f[c2] = V0 / t_rain * t[c2] # straight line within the duration of the rain, and 0 otherwise
    f[c3] = 0
    return f

What happens is that when I set the rain duration (t_rain) in dtdh2 to some large number, >9500, the solver converges and gives me the time required to empty the retention. However, when the time is smaller than that, is starts having problems at the beginning (near initial elevations), where the elevation (“time”) step becomes incredibly small, and the elevation values start to look like this: 234.99999996789282. I found this out by including a print(h) in dtdh2. The solver doesn’t terminate, instead it continues solving with the small “time” step, and thus never finishes (I interrupt it manually after some time). What can cause this problem, and what could I do to fix it?