# Need help with solve a system of delay differential equations

I am trying to solve and plot delayed SIRD model. My code is:

S4 =
NDSolve[
{x'[t] == - (0.3/(80*10^6)) x[t] y[t],
y'[t] == (0.3/(80*10^6)) x[t] y[t] -
(0.3/(80*10^6)) x[t - 14] y[t - 14] - 0.01 y[t],
z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - 14], a'[t] == 0.01 y[t],
x[0] == 80*10^6, y[t /; t <= 0] == E^t, y[0] == 150, z[0] == 0,
a[0] == 0 }, {x, y, z, a}, {t, 0, 200}]


In the solution I am getting the error message:

NDSolve::ndsz: At t == 164.8072392646089, step size is effectively zero; singularity or stiff system suspected.

And this results in values of graph that shouldnt be possible, active cases can't be negative.

What am I doing wrong?

I think there's a subtle inconsistency in your equations. You have to discount the recovery rate by the survival probability of infecteds, otherwise you double-count removals from y (both due to death and recovery), hence the negative y[t]. Luckily that is easy to do because survivorship is an exponential of the death rate and the infection time:

S4 = NDSolve[{
x'[t] == -(0.3/(80*10^6)) x[t] y[t],
y'[t] == (0.3/(80*10^6)) x[t] y[t] - (0.3/(80*10^6)) x[t - 14] y[t - 14] E^(-0.01*14) - 0.01 y[t],
z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - 14] E^(-0.01*14),
a'[t] == 0.01 y[t], x[0] == 80*10^6, y[t /; t <= 0] == E^t,
y[0] == 150, z[0] == 0, a[0] == 0}, {x, y, z, a}, {t, 0, 200}];

Plot[Evaluate[{x[t], y[t], z[t], a[t]} /. S4], {t, 0, 200}]


• ah thanks.it seems was indeed double counting recoveries. Nov 30, 2020 at 17:11

Remove y[0] == 150and set $$\tau=2$$ and replace $$y(t \cdots)$$ by $$y(t-\tau \cdots)$$ can remove some warm message.

τ = 2;
S4 =
NDSolve[{x'[t] == -(0.3/(80*10^6)) x[t] y[t - τ],
y'[t] == (0.3/(80*10^6)) x[t] y[t - τ] - (0.3/(80*10^6)) x[
t - 14] y[t - τ - 14] - 0.01 y[t - τ],
z'[t] == (0.3/(80*10^6)) x[t - 14] y[t - τ - 14],
a'[t] == 0.01 y[t - τ], x[0] == 80*10^6,
y[t /; t <= 0] == E^t, z[0] == 0, a[0] == 0}, {x, y, z, a}, {t, 0,
200}]
Plot[{x[t], y[t], z[t], a[t]} /. S4 // Evaluate, {t, 0, 200},
AxesOrigin -> {0, 0}, PlotRange -> All, AspectRatio -> 1]

• It still blows up, y` is gets very negative... (Version 12.1.1 Mac OS X.) Nov 30, 2020 at 14:06