How do you plot the solution of a system of odes against a parameter of the function and not time?

Hello,

I am trying to recreate some figures from a mathematical modelling paper, however I have arrived at one figure and unsure on how to plot it. This is the code that I have done for the function, it is an SIS model made up of two susceptible compartments.

function f=f(t,y)

f=zeros(3,1);

m=0.8;

A=200;

theta=0.004;

%lambda=0.00002; %R0>1

%lambda=0.00000989; %R0<1

n=0.7;

xi=0.012;

mu=0.035;

delta=0.3;

lambda1=0.1;

phi=0.01;

f(1)=m*A-theta*y(1)-lambda*y(1).*y(3)+n*xi*y(3)-mu*y(1);

f(2)=(1-m)*A+theta*y(1)-lambda*(1+delta*lambda1)*y(2).*y(3)+(1-n)*xi*y(3)-mu*y(2);

f(3)=lambda*y(1).*y(3)+lambda*(1+delta*lambda1)*y(2).*y(3)-(xi+phi+mu)*y(3);

This is the figure that I am wanting to replicate:

My main obstacle is trying to plot delta against the solution of y(:3), since it is not time. I have considered using a method, such as implict euler, however I am unsure on how to specify which is the argument of a function when the Matlab function has more inputs. I have tried using ode45 with this code

[t,u]=ode45(@f,[0 3000], [500 300 10]);

plot(t,u(:,3),'LineWidth',2)

but seem to have no luck in trying to change t to delta. Any help would be greatly appreciated.

Answer 1

[Cris LaPierre]

In your code, delta and m are constants. It is likely they ran multiple simulations to create that figure, varying the contants and capturing the desired value of y.

m = 0.2:0.2:0.8;
delta = 0:.05:2;
for a = 1:length(m)
    for b = 1:length(delta)
        [t,u]=ode45(@f,[0 3000], [500 300 10],[],m(a),delta(b));
        val(b,a) = u(end,3);
    end
end
plot(delta,val)
legend("m="+m','Location','Northwest')

function f=f(t,y,m,delta)
f=zeros(3,1);
% m=0.8;
A=200;
theta=0.004;
lambda=0.00002; %R0>1
%lambda=0.00000989; %R0<1
n=0.7;
xi=0.012;
mu=0.035;
% delta=0.3;
lambda1=0.1;
phi=0.01;
f(1)=m*A-theta*y(1)-lambda*y(1).*y(3)+n*xi*y(3)-mu*y(1);
f(2)=(1-m)*A+theta*y(1)-lambda*(1+delta*lambda1)*y(2).*y(3)+(1-n)*xi*y(3)-mu*y(2);
f(3)=lambda*y(1).*y(3)+lambda*(1+delta*lambda1)*y(2).*y(3)-(xi+phi+mu)*y(3);
end

[Jessica Rachel Furber]

This is great, thank you!