How to approximate impulsive differential equations?

I have a system of differential equations, one of which includes an impulse function. The system is seen below:

I cannot find any useful information on how to approximate this using ode45 (or other Matlab solvers) given the Dirac delta function in the problem.

The best idea I had was to try to break up my timespan into units of τ time blocks and then try to run ode45 on each time block, adding the b amount to B between each block. This looked closer to the right result, but still appears to fail.

For reference, if helpful, I am trying to replicate the results from the paper "Mathematical Model of Pulsed Immunotherapy for Superficial Bladder Cancer" from 2008, which can be found here: https://link.springer.com/article/10.1007/s11538-008-9344-z. I would appreciate any help anybody can give. Thank you.

[Paul]

This approach sounds correct in principal: break up my timespan into units of τ time blocks and then try to run ode45 on each time block, adding the b amount to B between each block.

Note that the you'll have to account for the increment of b at t = 0, i.e., the inittial conditions for the states are defined at t = 0-, and then B is incremented by b to get B(0+) (the other states don't change from t = 0- to t = 0+).

Will likely get more help by posting the code.

[Sam Chak]

Hi @Fiona

Are you aiming to generate an impulse train signal for the BCG instillation, similar to the one shown below? If that's the case, the question now is how to overcome the adaptive timestep of ode45 and instruct the solver to perform the instillation precisely every 7 days. I'm curious whether ode45 might miss the instillation at 6.999 days or 7.001 days.

By the way, you can find some good methods by @Paul to simulate the impulse response using ode45 in this thread.

tau = 7;                    % perform instillation every 7 days
f   = 1/tau;                % frequency of the instillation per day
fs  = f*1000;               % sample frequency
tEnd= 70;                   % end of simulation day
t   = 0:1/fs:tEnd;          % time span
Ins = zeros(size(t));       % initialization of Instillation
b   = 1;                    % dose injected
Ins(1:round(fs/f):end) = b;
plot(t, Ins), grid on, ylim([-0.5, 1.5])
xlabel Days, ylabel Dose
title('Instillation of  BCG')