Solving system of ODEs and Plotting

I am trying to graph the solution of a system of ODEs, given by

function system_of_DEs = f(t,x)
p = 0.5; A= 0.7; B= 0.75; C = 0.8;
h = x(1); l = x(2); v = x(3);
dhdt = B*(1-h-l)*h-v*h;
dldt = B*(1-h-l)*l + p*v*h-A*l;
dvdt = A*l+(1-p)*v*h-C*v;
system_of_DEs = [dhdt dldt dvdt]';
end

I want to plot the solution for C<1 for the time range 0-30 and then plot the solution to the same system of DE's on the same plot for C>1 for the time range 30-60.

Is there a way I can do this with one function? (Not have to create two seperate functions for the system of ODEs for C>1 and C<1). I've tried to do that here, but I got an error and I'm not sure why. Here's my function for the C>1 scenario.

function system_of_DEs_one = f(t,x)
p = 0.5; A= 0.7; B= 0.75; C = 1.1;
h = x(1); l = x(2); v = x(3);
dhdt = B*(1-h-l)*h-v*h;
dldt = B*(1-h-l)*l + p*v*h-A*l;
dvdt = A*l+(1-p)*v*h-C*v;
system_of_DEs_one = [dhdt dldt dvdt]';
end

Here's the .m script that I used to solve and plot these systems

timerange= [0 30];
timerange_1 = [30 60];
[t1, hlv1] = ode45(@system_of_DEs, timerange, [0.9,0.1,0.2]);
[t2, hlv2] = ode45(@system_of_DEs_one, timerange_1,hlv1(end,:));
plot(t1, hlv1(:,1), 'LineWidth', 2)
hold on
plot(t2, hlv2(:,1), 'LineWidth', 2)
hold off

I get an error running this and I'm not sure how to debug this.

I would also like if someone could find a way for me to only use one function/define one function for both scenarios described above, like a function that will take C as an input as well- f(t,x,C) and then call it. (I tried doing this and it didn't work--[t1, hlv1] = ode45(@system_of_DEs(1.1), timerange, [0.9,0.1,0.2]); it doesn't work (I put C=1.1, as an example).

Answer 1

[Alan Stevens]

Something like this?

timerange= [0 60];

[t, hvl] = ode45(@f, timerange, [0.9,0.1,0.2]);
plot(t, hvl, 'LineWidth', 2),grid
xlabel('t'),ylabel('h,v,l')
legend('h','v','l')

function system_of_DEs = f(t,x)
p = 0.5; A= 0.7; B= 0.75; C = [0.8 1.1];  %%%% C changed to a vector of values

ptr = 1; if t>30, ptr = 2; end

h = x(1); l = x(2); v = x(3);
dhdt = B*(1-h-l)*h-v*h;
dldt = B*(1-h-l)*l + p*v*h-A*l;
dvdt = A*l+(1-p)*v*h-C(ptr)*v;     %%%% C changed to C(ptr)
system_of_DEs = [dhdt dldt dvdt]';
end

Note that one should really use Event detection rather than using the "if" statement within the function. The results look ok to me, but only you can decide if they are acceptable.

[Maria Raheb]

This looks great- how can make the plot from 0-30 a different colour than the one for 30-60?

[Maria Raheb]

also I'm guessing that it's not possible to pass arguments into the function the ode45 solves, like ode45(@system_of_DEs(1.1), timerange, [0.9,0.1,0.2]).

[Alan Stevens]

Yes, it is possible to pass extra arguments to the ode function - see below, for example. Try

doc ode45

for more details.

The following also shows one way of plotting the two parts of the curve in different colours.

C = [0.8, 1];
timerange= [0 60];

[t, hvl] = ode45(@(t,x) f(t,x,C), timerange, [0.9,0.1,0.2]);

indx1 = find(t<=30); indx2 = find(t>=30);
t1 = t(indx1); t2 = t(indx2);
h1 = hvl(indx1,1); h2 = hvl(indx2,1);  % Similarly for v and l
plot(t1, h1,'r', t2, h2,'b'),grid
xlabel('t'),ylabel('h')

function system_of_DEs = f(t,x,C)
p = 0.5; A= 0.7; B= 0.75; 

ptr = 1; if t>30, ptr = 2; end

h = x(1); l = x(2); v = x(3);
dhdt = B*(1-h-l)*h-v*h;
dldt = B*(1-h-l)*l + p*v*h-A*l;
dvdt = A*l+(1-p)*v*h-C(ptr)*v;     %%%% C changed to C(ptr)
system_of_DEs = [dhdt dldt dvdt]';
end