Attempted to access betaI(2.72775e+06); index out of bounds because numel(betaI)=1

I have been working on this for a while now, but can't seem to figure out what's the problem. I keep getting the error:

    Attempted to access betaI(2.72775e+06); index out of bounds because numel(betaI)=1.

for the following code in MATLAB: (Function file followed by script)

    %System of differential equations

    function ydot = D12Eqns(t,y)

    global Lambda betaI betaA tau mu deltaIL nuIL muA deltaAL muT r2 r1 nuIS deltaSL nuIT h1 h2 a1 a2 

    XSS = y(1);
    XSL = y(2);
    XST = y(3);
    XIS = y(4);
    XIL = y(5);
    XIT = y(6);
    XAS = y(7);
    XAL = y(8);
    XAT = y(9);
    TSC = y(10);
    TIC = y(11);
    TAC = y(12);

    N = XSS+XSL+XST+XIS+XIL+XIT+XAS+XAL+XAT+TSC+TIC+TAC;

    ydot = zeros(12,1);

    ydot(1) = Lambda-betaI(XIS+XIL+XIT)*XSS/N-betaA(XAS+XAL+XAT)*XSS/N-tau(XST+XIT+XAT)*XSS/N-mu*XSS;
    ydot(2) = tau(XST+XIT+XAT)*XSS/N-betaI(XIS+XIL+XIT)*XSL/N-betaA(XAS+XAL+XAT)*XSL/N-(deltaSL+mu)*XSL-r1*XSL;
    ydot(3) = deltaSL*XSL-betaI(XIS+XIL+XIT)*XST/N-betaA(XAS+XAL+XAT)*XST/N-(mu+mu*T)*XST-r2*XST;
    ydot(4) = betaI(XIS+XIL+XIT)*XSS/N+betaA(XAS+XAL+XAT)*XSS/N-tau(XST+XIT+XAT)*XIS/N-(nuIS+mu)*XIS;
    ydot(5) = betaI(XIS+XIL+XIT)*XSL/N+betaA(XAS+XAL+XAT)*XSL/N+tau(XST+XIT+XAT)*XIS/N-(deltaIL+nuIL+mu)*XIL-h1*XSL;
    ydot(6) = betaI(XIS+XIL+XIT)*XST/N+betaA(XAS+XAL+XAT)*XST/N+deltaIL*XIL-(nuIT+mu+muT)*XIT-h2*XIT;
    ydot(7) = nuIS*XIS-tau(XST+XIT+XAT)*XAS/N-(mu+muA)*XAS;
    ydot(8) = nuIL*XIL+tau(XST+XIT+XAT)*XAS/N-(deltaAL+mu+muA)*XAL-a1*XAL;
    ydot(9) = nuIT*XIT+deltaAL*XAL-(mu+muA+muT)*XAT-a2*XAT;
    ydot(10)= r1*XSL+r2*XST-mu*TSC-betaI(XIS+XIL+XIT)*TSC/N-betaA(XAS+XAL+XAT)*TSC/N;
    ydot(11)= h1*XSL+h2*XIT-nuIS*TIC-mu*TIC;
    ydot(12)= a1*XAL+nuIS*TIC+a2*XAT-(muA+muT)*TAC;

    %Project Model Solver
    clear all

    global Lambda betaI betaA tau mu deltaIL nuIL muA deltaAL muT r2 r1 nuIS deltaSL nuIT h1 h2 a1 a2 

    Lambda = 6336000; %0.022 x N = 0.022 x 288000000%
    betaI = 0.170;
    betaA = 0.204;
    tau = 4;
    mu = 0.016;
    deltaIL = 0.05;
    muA = 0.5;
    muT = 0.2;
    deltaAL = 0.1;
    deltaSL = 0.0038;
    nuIS = 0.085;
    nuIT = 1;
    nuIL = 1.7;
    r1 = 3;
    r2 = 1;
    h1 = 3;
    h2 = 12;
    a1 = 12;
    a2 = 12;

    %%%== Initial conditions ==
    XSS0 = 168558250;
    XSL0 = 112668000; %XSL0-50000
    XST0 = 1968750; %XST0-50000
    XIS0 = 1587563;
    XIL0 = 1085500; %XIL0-25000
    XIT0 = 54687; %XIT0-25000
    XAS0 = 529188; %XAS0-12500
    XAL0 = 358000; %XAL0-12500
    XAT0 = 14062;
    TSC0 = 100000; %In TSC0, TIC0, and TAC0, we should reduce the numbers in the rest of the compartments (lines 27 to 35) and subtract from the total population and then divide them among the three compartments.
    TIC0 = 50000;
    TAC0 = 25000;

    y0 = [XSS0 XSL0 XST0 XIS0 XIL0 XIT0 XAS0 XAL0 XAT0 TSC0 TIC0 TAC0];

    tf = 40; %years;

    [t, y] = ode15s(@D12Eqns,[0 tf],y0);

    XSS = y(:,1);
    XSL = y(:,2);
    XST = y(:,3);
    XIS = y(:,4);
    XIL = y(:,5);
    XIT = y(:,6);
    XAS = y(:,7);
    XAL = y(:,8);
    XAT = y(:,9);
    TSC = y(:,10);
    TIC = y(:,11);
    TAC = y(:,12);

    plot(t,100.*(TSC)./N,'b')
    hold on
    axis square
    xlabel('Years')
    ylabel('Percentage of people in TSC out of whole population')

Would someone mind helping me? I'm working with a system of 12 ordinary differential equations. Thank you for your time and patience.

Answer 1

[Laurent]

Without running your code, I see that you define betaI as a scalar:

betaI = 0.170;

But then in your function you use betaI as if it is an array/vector:

ydot(1) = Lambda-betaI(XIS+XIL+XIT)*XSS/N-betaA(XAS+XAL+XAT)*XSS/N-tau(XST+XIT+XAT)*XSS/N-mu*XSS;

and so on.

You are trying to use the 'XIS+XIL+XIT'th value of betaI, while there is only 1, namely 0.170.

So either betaI is not correctly defined, or your ydot(1) and all the others that use betaI are not calculated correctly. Did you maybe forget a * between betaI and (XIS+XIL+XIT)?