Thanks, Daniel. And that is interesting. Your code produces correct results when the delta function is included in the single equation system. Mine apparently gets the sign wrong when the delta function is split out into a second equation.
I do get your point about numerically handling the delta function. The routine must understand when it has integrated over it, and I would think that would be troublesome for convergence. And it would of course be possible to use a delta function where it escaped from under the integral, and became a monster.
I also see that WhenEvent provides a perfectly good way of handling this. The only purpose of the delta under the integral is to increment the solution at that point, which can easily be done with WhenEvent. (WhenEvent was a great addition.)
For anyone interested, I ran into this just playing around. Suppose one takes a dose of medication every day at 10PM (2200). The dose moves to the bloodstream with a time constant of 1 hour. In the bloodstream, the concentration decays with a time constant of 7 hours. What is the quantity in the stomach and the blood over time?
eqs = {
(* stomach *)
s'[t] == -s[t]/ts,
(* bloodstream *)
b'[t] == s[t]/ts - b[t]/tb,
s[0] == 0,
b[0] == 0,
(* dose at 10PM *)
WhenEvent[Mod[t - 22, 24] == 0, s[t] -> s[t] + d]
};
{ss, bb} =
NDSolveValue[eqs /. {ts -> 1, tb -> 7, d -> 1}, {s, b}, {t, 0, 120}];
Plot[{ss[t], bb[t]}, {t, 0, 120}, PlotLegends -> {"s[t]", "b[t]"},
FrameLabel -> {"Time in hours", "Quantity or Concentration"},
PlotTheme -> "Scientific"]