Using a loop to get a script to repeat

Hi, I'm modelling disease spread using the SIR model and I've managed to code an m.file that works, but now I need to be able to get the model to loop, say 100 times, changing the value of the parameter beta by 0.01 each time to then get a range of values for R_0 which I can then put in a vector and plot to hopefully see a bifurcation.

Is there a way I can do this using loops rather than brute force going through each value of beta and then writing down the corresponding value of R_0?

This is my code:

function SIA

param.beta = 0.3;      % force of infection
param.mu = 0.5;        % birth rate
param.nu = 0;       % death rate
param.delta = 0.8;     % virulence of AIDS
param.gamma = 0.5;       % proportion who develop AIDS
param.alpha= 0;       % virulence of HIV

% Initial conditions are the values of all variables at time zero.

initial.S = 10^6;    % set the initial value of 'S'
initial.I = 10^5;     % set the initial value of 'I'
initial.A = 10^3;      % set the initial value of 'A'

inits = [initial.S; initial.I; initial.A]; % Initial conditions
end_time = 50;      % set how long to run for

function deriv = ode_system (t, x, param);

% Function to calculate derivatives of the SIR model

% Define N and calculate R_0

N = initial.S + initial.A + initial.I; % Population size
R_0 = param.beta / (param.gamma+param.nu+param.alpha)

% Define ODEs

S = x(1); % Number of susceptibles
I = x(2); % Number of HIV infected
A = x(3); % Number of AIDS infected

dS = +param.mu*N-param.beta * S * I/N-param.nu*S;
dI = +param.beta * S * I/N - (param.nu+param.gamma+param.alpha) * I;
dA = +param.gamma * I-(param.nu+param.delta)*A;

deriv = [dS; dI; dA];
end

% integrate the ODE system

[t, y] = ode45(@(t, x) ode_system(t, x, param),[0 end_time],inits,[]);

[Guillaume]

Please, remove all those blank lines that you probably spent a while putting in between each line of code and then click on the {} Code button to format them as code. Much faster and more readable this way.

Answer 1

[Tim]

Change "function SIA" to accept an input and return the integration. then call it from another script like so:

b=1:.01:whatever;
for i = 1:100
[t,y]=SIA(b(i));
result(i,1)= t;
result(i,2)= y;
end

and SIA will now be:

function [t,y] = SIA(b)

param.beta = b;
 blah
blah
blah