# NDSolve for solving delay equations with impulses

I want to solve a system of delay equations (DDE), that in some periods of time receive an impulse to one of its equations.

The equation (2) is the one receiving impulses of size $10^6$. So, I'm trying to use NDsolve of solve numerically this system with initial conditions $T=6*10^4$ ,$C_i=2*10^6$,$C_a=0$,$D=0$,$F_\beta = 0$ for $t<=0$.

dTdelay =
rt T[t] Log[K/T[t]] -
at Ca[t] T[t] (((1 - atb) etb)/(etb + Fb[t]) + atb);
ra Exp[-μci τ] Ci[t - τ] (Den[t - τ]/(
Den[t - τ] + θd)) +
re Exp[-μca τ] Den[t - τ] (Ca[t - τ]/(
Ca[t - τ] + θa)) - μca Ca[t];
dCidelay = -ra Exp[-μci τ] Ci[t - τ] (Den[t - τ]/(
Den[t - τ] + θd)) - μci Ci[t];
dFbdelay = rtb T[t] - Fb[t] μb;
dDendelay = (-Den[t] μd);

param3 = {rt -> 0.001060289, K -> 1*10^11, at -> 6*^-11, ht -> 5*10^8,
atb -> 0.69, etb -> 10^4,
re -> 6500, θd -> 10, θa -> 212,
ra -> 61, μca -> 0.01925, μci -> 0.007,
rtb -> 5.57*10^-6, μb -> 2.7, μd -> 0.009625};

dose = 10^6; τ = 265; tf = 1000;

teraphy2 = {WhenEvent[t == 168, Den[t] -> Den[t] + dose],
WhenEvent[t == 336, Den[t] -> Den[t] + dose],
WhenEvent[t == 504, Den[t] -> Den[t] + dose]};

sys = {T'[t] == dTdelay , Ca'[t] == dCadelay , Ci'[t] == dCidelay ,
Fb'[t] == dFbdelay , Den'[t] == dDendelay ,
T[t /; t <= 0] ==  6*^4, Ca[t /; t <= 0] ==  0,
Ci[t /; t <= 0] ==  2*^6, Fb[t /; t <= 0] ==  0,
Den[t /; t <= 0] ==  0};


Also i'm using WhenEvent to add the impulses at $t=168$ , $t=336$ , $t=504$.But NDsolve returns:

sol = NDSolve[
Join[sys, teraphy2] /. param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf}];

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


I know that the problems happen in $t=433$ when the delayed equations (3) and (4) start receiving the impulsed from $D[t-\tau]$ at $t=168$. Then suddenly (3) and (4) change. And that seems to cause problems to NDsolve.

NDsolve works fine until $t=433$ :

    sol = NDSolve[
Join[sys,
teraphy2, {WhenEvent[t == τ + 168, "StopIntegration"]} ] /.
param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf}];
Plot[Evaluate[{T[t], Ca[t], Ci[t], Fb[t], Den[t]} /. First[sol]], {t,
0, 433}, PlotRange -> All,
PlotLegends -> Placed[{"T", "Ca", "Ci", "Fβ", "D"}, Above],
Epilog -> {PointSize[0.02], Red, Apply[Point, {{168, 0}} ]} ]


Which returns:

I'm trying to overcome the problem integrating until $t=433$, save the data as interpolated functions. Using NDsolve to restart integrating at $t=433$ but im still receiving.

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


Any help of suggestion is welcome.

NDSolve will solve it if you decrease the AccuracyGoal but I am worried that your model may have some underlying problems. First, a solution up to t=500:

tf = 500;
sol = NDSolve[Join[sys, teraphy2] /. param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf},
AccuracyGoal -> 5];
Plot[Evaluate[{T[t], Ca[t], Ci[t], Fb[t], Den[t]} /. First[sol]], {t, 0, tf},
PlotRange -> All, PlotLegends -> Placed[{"T", "Ca", "Ci", "Fβ", "D"},
Above]]


That negative Ci worries me. Could be a numerical problem, but I think it might not be. Take a look at the first term in eqn 4, dCidelay. Ci decreases due to Ci τ time ago. Since Ci is decreasing with time, there is nothing preventing Ci from becoming negative, as it does. The following graphs support this argument:

Plot[Evaluate[{Ci[t], Ci[t - τ]} /. First[sol]], {t, 400, tf}, PlotRange -> {0, All}]
Plot[ra Exp[-μci τ] Ci[t - τ] (Den[t - τ]/(Den[t - τ] + θd)) /.
param3 /. sol[[1]], {t, 400, tf}]


If you let the next pulse hit, all hell breaks loose:

tf = 710;
sol = NDSolve[Join[sys, teraphy2] /. param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf}, AccuracyGoal -> 5];
Plot[Evaluate[{T[t], Ca[t], Ci[t], Fb[t], Den[t]} /. First[sol]], {t, 400, tf},
PlotRange -> All, PlotLegends -> Placed[{"T", "Ca", "Ci", "Fβ", "D"}, Above]]


So, my guess is it's more of a model problem than an NDSolve problem.