Varying parameters in ODE and plot the different outcomes

I have four differential equations:

And my parameters are:

h=2; tau=1.5; rho=0.07; v=0.002; u=0.8; R=0.1; r=R*(1-exp(-v*t.^2);
% I have set r = R = 0.1 and thus r = R*(1-exp(-v*t.^2)

I want to vary h, tau, rho, v and R individually for a few values below and above baseline value, about up to 50% while the remaining parameters not being tested are held constant. I then want to plot each test on a single plot along with the baseline, so I should have five plots in total.

I have never done anything like this (and tried to follow an example in the book "Undergraduate Texts in Mathematics" by Shonkwiler and Herod on page 387), and my approach was to use ode45 but I don't seem to be able to get it to work for a varying parameter. This is my code so far:

for i=1:3;
    i=i+1;
S0=1; I10=0; I20=0; R0=0;
h=[1 1.5 2 2.5 3]; tau=1.5; rho=0.07; v=0.002; u=0.8; 

dt=0.01;
t = 0:dt:50;
timeSpan = [0 50];

R=0.1; r=R*(1-exp(-v.*t.^2));
beta=h.*exp(-h.*tau)/(1-exp(-h.*tau));

odes= @(t,x) [-h.*x(1)+rho*x(2)+beta*x(4); 
      (1-u)*h.*x(1) - rho*x(2) - u*h.*x(2); 
      u*h.*x(1) + u*h.*x(2) - r(i)*x(3);
      r(i)*x(3) - beta*x(4)];
  
[tut xut]=ode45(odes,timeSpan,[S0 I10 I20 R0]);

plot(tut,xut)
end

I get the error message: Error using vertcat

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in
testtt>@(t,x)[-h.*x(1)+rho*x(2)+beta*x(4);(1-u)*h.*x(1)-rho*x(2)-u*h.*x(2);u*h.*x(1)+u*h.*x(2)-r(i)*x(3);r(i)*x(3)-beta*x(4)]
(line 16)
odes= @(t,x) [-h.*x(1)+rho*x(2)+beta*x(4);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options,
  varargin);

Error in testtt (line 21)
[tut xut]=ode45(odes,timeSpan,[S0 I10 I20 R0]);

Could this be solved by ode45 or should I do something else? All help is truly appreciated!

Answer 1

[Alan Stevens]

More like this perhaps:

timeSpan = [0 50];
%     S  I1 I2 R  S  I1 I2 R  S  I1 I2 R  
P0 = [1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0];
% 3 sets of initial values: one for hi, baseline and lo respectively

selectiontype = [' h '; 'tau'; 'rho'; ' v '; ' u '];
select = 1;  % 1 for h, 2 for tau, etc.

% call ode solver
[t, P]=ode45(@dPdtfn,timeSpan,P0,[],select);

% plot results
% 3 sets of 4 columns of P, for 'hi', 'baseline' and 'lo' respectively.
subplot(3,1,1)
plot(t,P(:,1:4)),grid
xlabel('time'),ylabel(['Population (hi ',selectiontype(select,:),')'])
legend('S','I1','I2','R')
subplot(3,1,2)
plot(t,P(:,5:8)),grid
xlabel('time'),ylabel(['Population (base ',selectiontype(select,:),')'])
legend('S','I1','I2','R')
subplot(3,1,3)
plot(t,P(:,9:12)),grid
xlabel('time'),ylabel(['Population (lo ',selectiontype(select,:),')'])
legend('S','I1','I2','R')


function dPdt = dPdtfn(t,P,select)

        % Baseline values
        hb=2; taub=1.5; rhob=0.07; vb=0.002; ub=0.8;
        
        % Default hi and lo values
        hhi = hb; tauhi = taub; rhohi = rhob; vhi = vb; uhi = ub;
        hlo = hb; taulo = taub; rholo = rhob; vlo = vb; ulo = ub;
        
        % Modify appropriate hi and lo values
        if select == 1
           hhi = 1.5*hb;  hlo = 0.5*hb;
        elseif select == 2
            tauhi = 1.5*taub; taulo = 0.5*taub;
        elseif select == 3
            rhohi = 1.5*rhob;  rholo = 0.5*rhob;
        elseif select == 4
            vhi = 1.5*vhi; vlo = 0.5*vb;
        else
            uhi = 1.5*ub; ulo = 0.5*ub;
        end
            
        betahi=hhi.*exp(-hhi.*tauhi)/(1-exp(-hhi.*tauhi));
        betab=hb.*exp(-hb.*taub)/(1-exp(-hb.*taub));
        betalo=hlo.*exp(-hlo.*taulo)/(1-exp(-hlo.*taulo));  
        rhi = 0.1*(1-exp(-vhi.*t.^2));
        rb = 0.1*(1-exp(-vb.*t.^2));
        rlo = 0.1*(1-exp(-vlo.*t.^2));
        
        % Extract population types from P 
        S = P(1); I1 = P(2); I2 = P(3); R = P(4);
        
        % Rates of change
        dPdthi = [-hhi.*S+rhohi*I1+betahi*R; 
                (1-uhi)*hhi.*S - rhohi*I1 - uhi*hhi.*I1; 
                uhi*hhi.*S + uhi*hhi.*I1 - rhi*I2;
                rhi*I2 - betahi*R];
        dPdtb = [-hb.*S+rhob*I1+betab*R; 
                (1-ub)*hb.*S - rhob*I1 - ub*hb.*I1; 
                ub*hb.*S + ub*hb.*I1 - rb*I2;
                rb*I2 - betab*R];
        dPdtlo = [-hlo.*S+rholo*I1+betalo*R; 
                (1-ulo)*hlo.*S - rholo*I1 - ulo*hlo.*I1; 
                ulo*hlo.*S + ulo*hlo.*I1 - rlo*I2;
                rlo*I2 - betalo*R];
            
            dPdt = [dPdthi;
                    dPdtb;
                    dPdtlo];
        
end