Heaviside function in ODE

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.

Answer 1

[Torsten]

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.

[Arjun Chakrawal]

Pascal Schulthess-- could you please share your code on how did you actually implemented the suggestions by Torsten?

[Sayan Batabyal]

@Torsten, Do you know how to change the value of flag inside the events function? If you can, please elaboarte with a small code snippet.