My code is used to simulate various forms of the differential equation -2 / (h(t) - λ(t)) for some starting point h0 and driver λ(t). I’ve ran this for some h0 and λ(t) with good success (such as λ(t) = 0 and λ(t) = sqrt(t)), however I noticed that in some situations (such as this current one) the solution can veer off quite a bit. In this case here, with λ(t) = t, the differential equation to solve is dhdt = -2 / (h(t) - t) with h0 = -2 .

def dh_dt_1cc(t, h, lambda1):
    return (-2) / (h - (lambda1(t)))

def solve_diffeq_1cc(h0, t, lambda1):
    sol = solve_ivp(dh_dt_1cc, t_span=[0, t], y0=[h0], args=(lambda1,), events=event, rtol=1e-5, atol=1e-9, max_step=1e-3)
    return sol

#event is used to stop calculating the solution when the solution intersects the lambda,
#the hitting-time between the solution and the lambda is the result of interest
def event(t, h, lambda1):
    return h[0] - lambda1(t)

event.terminal = True
event.direction = 0

if __name__ == "__main__":
    h0 = -2
    lambda1 = lambda t: t
    sol = solve_diffeq_1cc(harr, 100, lambda1)

The analytical solution provides h(t) = t - 2, which means that it will never intersect with the lambda. This is fine, however when I take a look at the solution generated by solve_ivp using the RK45 method, I notice that the solution stays mostly correct up until time 55 where the error between the numerical and the analytical solution slowly starts to increase. This makes me worry that the solution may be inaccurate when I run it for other drivers that may have and intersection between the solution and the driver. How can I prevent this error from arising? How can I quantify the error that the solution may generate when I am unable to compare it against an analytical solution? Below is a plot with error between the numerical solution and the analytical solution.

I’ve tried using the implicit method Radau to no avail.