DDE package is a fast library written in Haskell that I have extracted after my PhD. Previously, my work involved a lot dynamical systems described by DDEs. I kept switching between several tools, still I had no fast and flexible library. So here is one now.

Below is an example of a dual-delay model taken from a recent publication. It is given by the following delay differential equation:

$$ \varepsilon \frac{dx(t)}{dt}+x(t)+\delta\cdot\intop_{t_{0}}^{t}x(\xi)d\xi = (1-\gamma)f\left(x_{\tau_{1}}\right)+\gamma f\left(x_{\tau_{2}}\right), $$

where $ \tau_2 = 100 \tau_1$ are delay times and $\Phi_0$ is a control parameter. Let us fix dynamics parameters $\varepsilon = 0.01$, $\delta = 0.009$, and $\gamma = 0.5$.

We transform the integro-differential equation into an equivalent system of equations:

$$\begin{aligned}
\varepsilon \dot x &= -x -\delta\cdot y + (1-\gamma)f \left( x_{\tau_1} \right) + \gamma f \left( x_{\tau_2} \right),\\

\dot y &= x.
\end{aligned}
$$

Now, this system can be simply written in a vectorial (V2) form:

```
-- Equation right-hand side
rhs phi0 = DDE.RHS derivative
where
derivative ((V2 x y), (DDE.Hist histSnapshots), _) = V2 x' y'
where
-- DDE Eq. (4) from arXiv:1712.03283
x' = (-x - delta * y + (1 - gamma) * f x_tau1 + gamma * f x_tau2) / epsilon
y' = x
f = airy phi0
-- Delay terms where tau2 / tau1 = 100
(V2 x_tau1 _):(V2 x_tau2 _):_ = histSnapshots
-- Constants
epsilon = 0.01
gamma = 0.5
delta = 0.009
```

We define the nonlinear Airy function $f(x) = \beta [1 + m \sin^2(x + \Phi_0)]^{-1}$, $\beta = 1.6$, $m = 50$ as follows:

```
airy phi0 x = beta / (1 + m * (sin (x + phi0))^2)
where
m = 50
beta = 1.6
```

The complete model is available here.