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.
Best Answer