MATLAB: Solving a system of differential equations using ODE45 by switching between two functions depending upon the solution obtained in last iteration.

conditional statementdifferential equationsMATLAB

Problem: I have a system of two coupled differential equations but i want to put an conditional statement so that depending upon that condition it evaluates the system with different parameter values.
My attempt: I couldn't do the problem and hence I don't have a MWE but here is what i tried. I defined two separate functions with the desired parameters and somehow I wanted to check for the result i obtained in the last iteration and then using an if block to satisfy my condition i wanted to let the solver know which function to evaluate.
I am aware about these event functions and examples like ballode, but i don't know how to formulate it in that way.
timerange= 0:0.5:200;
IC= [0.1,0.1];%initial conditions
threshold=0.5;
% [t,y] =ode45(@(t,y) fn_1(t,y),timerange, IC);
%Not a MWE but more of how i want the system to behave.
%I want that whenever the solution of first DE exceed the threshold,
%It switches to the fn_2
if y(:,1)>=threshold
    [t,y] =ode45(@(t,y) fn_2(t,y),timerange, IC);
else
    [t,y] =ode45(@(t,y) fn_1(t,y),timerange, IC);
end
%---
    
plot(t,y(:,1),'r');
hold on
plot(t,y(:,2),'b');
xlabel('n')
ylabel('I')
grid on
function rk1 =fn_1(t,y)
r      = 0.5; K = 0.5;
eps    = 0.001;  rho= 0.02;
alpha1 = 2.15; c1= 0.025;
gammaI = 0.02; I0= 0.01;
n= y(1);
I= y(2);
    
rk1(1)= r*n*(1- n/K)-eps*n*I;
rk1(2) = I0 + (rho*I*n)/(alpha1+n) -c1*I*n - gammaI*I;
rk1=rk1(:);
end
function rk1 =fn_2(t,y)
r      = 0.5;   K= 0.5;
eps    = 0.001; rho= 0.02;
alpha1 = 2.15;  c1 = 0.025;
I0     = 0.5;     gammaI = 0.2;
n= y(1);
I= y(2);
    
rk1(1)= r*n*(1- n/K)-eps*n*I;
rk1(2) = I0 + (rho*I*n)/(alpha1+n) -c1*I*n - gammaI*I;
rk1=rk1(:);
end
Any help with the code or a hint in the right direction would be very helpful.
Thank you.

Best Answer

According to the problem you defined in your comment, this is the correct way to implement it using ode45
timerange = 0:0.5:200;
IC= [0.1,0.1];%initial conditions
threshold=0.5;
[t,y] =ode45(@(t,y) fn_1(t,y), timerange, IC);
plot(t,y(:,1),'r');
hold on
plot(t,y(:,2),'b');
xlabel('n')
ylabel('I')
grid on
function rk1 =fn_1(t, y)
n = y(1);
I = y(2);
if n < 0.1
    r      = 0.5; K = 0.5;
    eps    = 0.001;  rho= 0.02;
    alpha1 = 2.15; c1= 0.025;
    gammaI = 0.02; I0= 0.01;
else
    r      = 0.5;   K= 0.5;
    eps    = 0.001; rho= 0.02;
    alpha1 = 2.15;  c1 = 0.025;
    I0     = 0.5;     gammaI = 0.2;
end
rk1(1) = r*n*(1- n/K)-eps*n*I;
rk1(2) = I0 + (rho*I*n)/(alpha1+n) -c1*I*n - gammaI*I;
rk1 = rk1(:);
end
However, even the 2nd set of parameters does not decrease the value of n.
Related Question