how to make a for loop in an ODE 45

data = [0 86; 1 86; 2 117; 6 120; 7 130; 8 135; 13 169; 16 179; 23 224; 27 230; 29 242; 36 234; 41 244; 50 245; 59 270; 63 271; 64 309; 69 354; 72 438; 77 476; 78 508; 85 528; 91 599; 99 759; 104 779; 105 844; 111 888; 113 964; 118 1048; 121 1093; 125 1201; 128 1322; 131 1437; 132 1617; 136 1766; 140 1835; 141 1963; 143 2115; 147 2225; 149 2458; 150 2599; 156 3052; 165 3944; 167 4269; 171 4366; 175 4963; 177 5325; 181 5843; 183 6242; 185 6553; 190 7157; 192 7470; 197 8011; 199 8376; 204 8973; 206 9191; 211 9911;214 10114; 218 13676; 220 13540; 225 13015; 227 13241; 232 14068; 234 14383; 239 15113; 241 15319; 246 15901; 248 16899; 253 17111; 260 17908; 267 18569; 274 19463; 281 20171; 288 20712; 295 21261; 302 21689; 309 22057; 316 22460; 323 22859; 330 23218; 337 23694; 344 23934; 351 24247; 358 24666; 365 24872; 372 25178; 379 25515; 386 25791; 393 26044; 400 26277; 407 26593; 414 26724; 421 26933; 428 27013; 435 27145; 442 27237; 449 27305; 456 27443; 463 27514; 470 27573; 477 27642; 484 27705; 491 27748; 498 27862; 505 27929; 512 27952; 519 28005; 527 28073; 534 28147; 541 28220; 548 28295; 555 28388; 562 28421; 569 28454; 576 28476; 583 28539; 590 28571; 596 28599; 603 28598; 610 28601; 617 28601; 624 28601; 631 28604; 638 28601; 645 28601; 652 28601; 659 28601];

function Math462hw3Q5parta

alpha=0.05;

beta=0.2;

gamma=0.125;

N=50000;

S=1-I-E;

for i= 0:1:659

I=data(i,2)/N;

end

E=2*I;

R=0;

y=N*alpha*E;

vector=[S,E,I,R,y];

tspan=0:1:659;

[t,soln]=ode45(@(t,vector)dotfunctions(t,vector,alpha,beta,gamma),tspan,vector);

%plot function

plot(t,soln);

xlabel('time');

ylabel('value');

title('SIERy model');

function ddt=dotfunctions(~,siery,alpha,beta,gamma)

S=siery(1);

I=siery(2);

R=siery(3);

E=siery(4);

y=siery(5);

Sdot=-beta*S*I;

Edot=beta*S*I-alpha*E;

Idot=alpha*E-gamma*I;

Rdot=gamma*I;

ydot=N*alpha*E;

ddt=[Sdot;Idot;Rdot;Edot;ydot];

end

end

This is my code for simulating a model using ODE45, here are two questions:

1. Can I write a for loop in the first few lines just to loop through all values of i
2. Just wondering can someone help me figure out why this code is not running.

Answer 1

[Alan Stevens]

More like this I think:

data = [0 86; 1 86; 2 117; 6 120; 7 130; 8 135; 13 169; 16 179; 23 224; ...
    27 230; 29 242; 36 234; 41 244; 50 245; 59 270; 63 271; 64 309; ...
    69 354; 72 438; 77 476; 78 508; 85 528; 91 599; 99 759; 104 779; ...
    105 844; 111 888; 113 964; 118 1048; 121 1093; 125 1201; 128 1322; ...
    131 1437; 132 1617; 136 1766; 140 1835; 141 1963; 143 2115; 147 2225; ...
    149 2458; 150 2599; 156 3052; 165 3944; 167 4269; 171 4366; 175 4963; ...
    177 5325; 181 5843; 183 6242; 185 6553; 190 7157; 192 7470; 197 8011; ...
    199 8376; 204 8973; 206 9191; 211 9911;214 10114; 218 13676; 220 13540;...
    225 13015; 227 13241; 232 14068; 234 14383; 239 15113; 241 15319; ...
    246 15901; 248 16899; 253 17111; 260 17908; 267 18569; 274 19463; ...
    281 20171; 288 20712; 295 21261; 302 21689; 309 22057; 316 22460; ...
    323 22859; 330 23218; 337 23694; 344 23934; 351 24247; 358 24666; ...
    365 24872; 372 25178; 379 25515; 386 25791; 393 26044; 400 26277; ...
    407 26593; 414 26724; 421 26933; 428 27013; 435 27145; 442 27237; ...
    449 27305; 456 27443; 463 27514; 470 27573; 477 27642; 484 27705; ...
    491 27748; 498 27862; 505 27929; 512 27952; 519 28005; 527 28073; ...
    534 28147; 541 28220; 548 28295; 555 28388; 562 28421; 569 28454; ...
    576 28476; 583 28539; 590 28571; 596 28599; 603 28598; 610 28601; ...
    617 28601; 624 28601; 631 28604; 638 28601; 645 28601; 652 28601; 659 28601];
[nrows, ncols] = size(data);

alpha=0.05;
beta=0.2;
gamma=0.125;
N=50000;
for i= 1:nrows  % indices start with 1 in Matlab, not 0
    I=data(i,2)/N;
end
E=2*I;
S=1-I-E;  % needs to come after I and E are defined
R=0;
y=N*alpha*E;
vector=[S,E,I,R,y];
tspan=0:1:659;
[t,soln]=ode45(@(t,vector)dotfunctions(t,vector,alpha,beta,gamma,N),tspan,vector);
%plot function
subplot(2,1,1)       % separate plots as y scale swamps the others
plot(t,soln(:,1:4)), grid
xlabel('time');
ylabel('value');
legend('S','E','I','R')
title('SIERy model');
subplot(2,1,2)
plot(t,soln(:,5)), grid
xlabel('time');
ylabel('y value');


    function ddt=dotfunctions(~,siery,alpha,beta,gamma,N)  % need to pass in N as well
        S=siery(1);
        I=siery(2);
        % R=siery(3);
        E=siery(4);
        % y=siery(5);
        
        Sdot=-beta*S*I;
        Edot=beta*S*I-alpha*E;
        Idot=alpha*E-gamma*I;
        Rdot=gamma*I;
        ydot=N*alpha*E;
        ddt=[Sdot;Idot;Rdot;Edot;ydot];
    end