I made the following MATLAB code using symbolic toolbox and it represents drug dosage dissolving in the body. Now I want to consider a patient who takes multiple doses of the same capsule taken at t=0, t-120, and t=240 min. How would I go about updating my symbolic solver to get this new output. I know that I have to work with the dirac(t) function but I am not sure how to go about implementing it. Any ideas? Thanks in advanced.

MATLAB CODE

%clear screen

clc;

%use "syms" command to identify specific variables and functions as

%symbolic

syms k_a k_d k_m D_0 r;

% use "syms" command to also specify functions

syms B(t) D(t);

%Input equation

f(t) = D_0*(1-exp(-r*t));

%find derivatives

D1B = diff (B, t, 1); % first derivative of B(t)

D1D = diff (D, t, 1); % first derivative of D(t)

%g - forcing function

g = diff (f, t, 1); % first derivative of the input function

%ODEs combined with "odes" command

odes = [odeD; odeB];

%Initial conditions

C1 = D(0) == 0;

C2 = D(0) == 0;

%Use "dsolve" command to solve ODEs symbollically

y = dsolve (odes,[cond1;cond2]);

%Assign Numeric Values

k_a = 0.025; % (min^-1 - absorption rate)

k_d = 0.005; % (min^-1 - digestion rate)

k_m = 0.015; % (min^-1 - metabolic rate)

D_0 = 50; % mg - initial dose

r = 0.1151; %(min^-1) - Drug release rate

% Define plotting range

t = linspace(0,6*60,2000);

% Use "eval" command for "yNum"

Dnum = eval(y.D);

Bnum = eval(y.B);

W = eval(f) - (Bnum + Dnum);

%Plot

plot(t,Dnum,t,Bnum,t,W);