I have an inseparable, first-order ODE of the form $$f'(t)=g(f'(t), f(t))$$
where $g$ is a rather monstrous combination of terms making it impossible to find some $h$ such that $f'(t)=h(f(t))$. I'd like to numerically integrate $f'(t)$ over time, ideally with something like scipy's solve_ivp(). However, I'm not sure how to best approach this problem as solve_ivp() (and the other ODE solvers I'm aware of) require the input function to have the form $f'(t)=h(f(t))$.
An option that comes to mind is to brute force $f(0)$ and $f(\Delta t)$ to estimate $f'(0)=\frac{f(\Delta t)-f(0)}{\Delta t}$. With $f'(0)$ estimated, I could then pass that into solve_ivp() and update this $f'(t)$-estimate each time solve_ivp() recalculates the derivative, passing the updated version to the subsequent call. I don't expect $f(t)$ to change rapidly with $t$, so using a "stale" value for $f'(t)$ is probably ok, but it feels like an awkward/wrong way of doing things.
Is there a simpler/better approach I've overlooked? Are there other Python packages better suited to this type of problem?