MATLAB: Heaviside function in ODE

MATLABode

I have the following set of ODEs
dxdt(1) = -ka*x(1);
dxdt(2) = ka/V*x(1) - ke*x(2);
dxdt(3) = -a1*exp(-b1*x(3)).*x(3) + (x(2) - c1)*H(x(2) - c1);
wherein H(x) = 1 for x >= 0 and 0 otherwise is a Heaviside function. I've learned that such an ODE should be solved by using events.
This is what I currently have:
function mcm
% Time
tstart = 0;
tfinal = 2*24;
tspan = tstart:0.1:tfinal;
% Initial conditions
x0 = [200; 0; 0];
% Parameters
ka = 2.4;
V = 14;
ke = 0.39;
a1 = 0.7;
b1 = 0.31;
c1 = 3.7;
% Options with event function
options = odeset('Events', @heavi);
% Solve until first terminal event
[t, x, te, xe, ie] = ode45(@model, tspan, x0, options);
% Plot
plotnames = ['x1'; 'x2'; 'x3'];
for i = 1:length(x0)
    subplot(3, 1, i);
    plot(t, x(:, i))
    title(plotnames(i, :));
end
% Model
    function dxdt = model(t, x)
          % Initialise
          dxdt = zeros(3, 1);
          % Model equations
          dxdt(1) = -ka*x(1);
          dxdt(2) = ka/V*x(1) - ke*x(2);
          dxdt(3) = -a1*exp(-b1*x(3)).*x(3);
      end
% Event function
    function [value, isterminal, direction] = heavi(t, x)
        value = x(2) - c1;
        isterminal = 0;
        direction = 0;
    end
end
But I'm not sure how to add the Heaviside part to the third ODE. Any help would be greatly appreciated.

Best Answer

1. Set isterminal=1 in the event function and organize your code as in the "ballode" example.
2. Switch a variable "flag" in "mcm" from 0 to 1 and vice versa each time you return from ode45 (i.e. when the event function detects a zero crossing of "value") and set
dxdt(3) = -a1*exp(-b1*x(3)).*x(3)+flag*(x(2)-c1);
Best wishes
Torsten.
Related Question