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

8 views (last 30 days)

Show older comments

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.

##### 7 Comments

Ameer Hamza
on 4 May 2020

### Accepted Answer

Ameer Hamza
on 4 May 2020

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.

##### 4 Comments

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!