# Numerically solve delay ODE

I am trying to numerically solve a simple delay ODE in mma for the first time. I think this is likely a simple fix.

I first tried to solve system as ordinary differential equation system for up to 𝑡=l. And then use these solutions as history functions for the delay differential equation system. I keep getting the error: NDSolve::dvnoarg The function na appears with no arguments

ClearAll["Global*"]
c = 200; Q = 4; P = 1.92*10^-6; ea = 2.6*10^-6; ga =
6.24*10^-8; l = 0.5; ba = 98; pr = 0.2;
sol0 = NDSolve[{
r'[t] == pr (c - r[t]) - P r[t]/(Q + r[t]) (na[t] + ma[t]),
na'[t] == P r[t]/(ea (Q + r[t])) na[t] - pr na[t] - ga
na[t] p[t],
ma'[t] == ga na[t] p[t] - pr ma[t],
p'[t] == -pr p[t] - p[t] ga na[t] ,
r == 1000, ma == 0, na == 100, p == 1000},
{r, na, ma, p}, {t, 0, l}][];

sol = NDSolve[{
r'[t] == pr (c - r[t]) - P r[t]/(Q + r[t]) (na[t] + ma[t]),
na'[t] == P r[t]/(ea (Q + r[t])) na[t] - pr na[t] - ga
na[t] p[t],
ma'[t] == ga na[t] p[t] - pr ma[t] -
Exp[-pr l] ga na[t - l] p[t - l],
p'[t] == ba Exp[-pr l] ga na[t - l] p[t - l] - pr p[t] -
p[t] ga na[t] ,
r == 1000, na[t /; t <= l] == sol0[], ma == 0,
p[t /; t <= l] == sol0[]},
{r, na, ma, p}, {t, 0, 100}][];


Any suggestions are much appreciated.

• You might look at sol0[] and sol0[]. You probably want (na[t] /. sol0), or to use NDSolveValue for sol0 (the first NDSolve call). Apr 14, 2019 at 15:12
• Thanks, (na[t] /. sol0) worked. Apr 14, 2019 at 23:22
• BTW, what's the model? Looks like a nutrient-consumer-virus system. Apr 15, 2019 at 7:39
• You are correct @Chris K Apr 15, 2019 at 14:26
• Nice, pretty close to my own research. I added an "ecology" tag. Good luck. Apr 15, 2019 at 14:44

I'm not sure the separate NDSolve to initialize the initial conditions is even necessary. If you just use your second NDSolve with the initial conditions of the first, it gives a warning NDSolve::ihist but seems to run:

sol = NDSolve[{
r'[t] == pr (c - r[t]) - P r[t]/(Q + r[t]) (na[t] + ma[t]),
na'[t] == P r[t]/(ea (Q + r[t])) na[t] - pr na[t] - ga na[t] p[t],
ma'[t] == ga na[t] p[t] - pr ma[t] - Exp[-pr l] ga na[t - l] p[t - l],
p'[t] ==  ba Exp[-pr l] ga na[t - l] p[t - l] - pr p[t] - p[t] ga na[t],
r == 1000, na == 100, ma == 0, p == 1000},
{r, na, ma, p}, {t, 0, 400}, AccuracyGoal -> 128][];

LogPlot[Evaluate[{na[t], ma[t], p[t]} /. sol], {t, 0, 400},
PlotRange -> All] However note that I had to really crank up AccuracyGoal due to the very low levels that na[t]` reaches. Maybe other parameter values would be less violent.