[Feature request] A delay differential equations solver

Dear Patrick,

Thanks for this library, it's pretty neat!

I'm currently supervising a bachelor's students thesis on NeuralODEs, and we've been meaning to use diffrax to study delay differential equations. I'm currently trying to implement a delay differential equations solver inside diffrax, and I would like to know if I'm in the right track.

For context, let me give a brief overview of how delay differential equations work, and how such a solver could be implemented. In it's simplest form, a (constant) Delay Differential Equation (DDE) has a vector field $f$ that depends not only on the current state $y(t)$, but also on $y(t-\tau)$ where $\tau\in\mathbb{R}_{>0}$. In other words

$$y'(t) = f(t, y(t), y(t-\tau)).$$

Initial value problems involving DDEs provide a history instead of a single initial value ($y(t) = \phi(t)$ for $t \leq 0$, for example), and are solved in chunks using the "method of steps". Shortly put, one solves an IVP in intervals of the form $[t_0 + k\tau, t_0 + (k+1)\tau]$. (More details in Chap. 9 of this reference).

In practice, solving a DDE numerically can be done by selecting the right step-size such that $y(t-\tau)$ is always in the grid. To predict y_{t+1} we need to evaluate the vector field in y_{t} and in some y_{t-k} corresponding to $y(t-\tau)$. Other ways of doing it would involve e.g. building an Hermite interpolation between the two relevant points in the grid if $y(t-\tau)$ happens to lie outside of the grid, but I plan to focus on the first alternative.

How could I adapt diffrax to let me pass terms of the form f(t, y(t), y(t-tau))dt? I imagine I have to implement a new DelayTerm that inherits from AbstractTerm with a different vf method; since solvers need to evaluate these vector fields, I imagine I would also need to modify/create a new one in which vf is called with the right signature, right?

I'm of course happy to contribute the implementation to diffrax once it's up and running.

Start by reading the AbstractTerm and solver interfaces, then examine how vf is passed and evaluated. The proposal mentions a new DelayTerm and solver changes, but no files or tests are identified. Done would mean a defined DDE API, solver behavior for delayed state access, and tests covering the supported delay-solving approach.