It looks like you have a simple clearance model, and at each time in a
, you want the dose to be delivered instantaneously. That is, at each time in a
, the amount of drug in the subject increases by 57.947.
If that is the model that you have in mind, implementing it in the formula for xdot
will not work well. You would actually need to implement it as a "delta function", and lsode
won't work with that.
Instead, I recommend solving the problem in stages, corresponding to the time intervals [0, 1], [1, 7], [7, 21], and [21, 28]. During one stage, the differental equation is simply xdot = -0.4077*x
. In the first stage, the initial condition is 0. In the next stage, the initial condition is the final value of the previous stage plus the dosage amount 57.947.
Here's a script:
dose = 57.947;
a = [1 7 21 28 30];
x0 = 0.0;
t0 = 0;
t = [];
y = [];
for t1 = a
tinterval = linspace(t0, t1, 2*(t1 - t0) + 1);
yinterval = lsode(@(x, t) -0.4077*x, x0, tinterval);
t = [t tinterval];
y = [y yinterval'];
t0 = t1;
x0 = yinterval(end) + dose;
end
plot(t, y, "linewidth", 2);
The script creates this plot:
Note that the differential equation xdot = -k*x
has the solution x0*exp(-k*(t-t0))
, so the call to lsode
could be replaced with
yinterval = x0*exp(-0.4077*(tinterval - t0));
If you do this, also remove the transpose from yinterval
two lines below:
y = [y yinterval];
If you want to keep the implementation of the administration of the drug dose within the formula for xdot
, you'll need to distribute it over a small time interval. It could be implemented as a short rectangular pulse of width w
and height 57.974/w
. You'll also need to ensure that lsode
takes sufficiently small internal time steps (less than w
) so that it "sees" the drug dose.