MATLAB: What can I do to avoid the error “Index exceeds matrix dimensions ” in the simulation for systems of equation

numerical simulation of ode system

    mu = pars(1);            pi = pars(2);          omegaA = pars(3);

Best Answer

Hi,
gammacAc=0.1128;
was missing in parameter declaration and in pars. Try:
function Meningitis_Coinfection
    mu=1/(56*365);
    pi=mu*182202000*(1/365);
    omegaA=(1/5)*(1/365);
    kappaA = 0.0031;
    kappac = 0.0031;
    kappaAc = 0.0031;
    nuA=0.14868; 
    epsilon=((1.90)/2);
    sigmaA=0.0548;
    sigmac=0.0548;
    sigmaAc=0.0548;
    gammacA= 0.1118;
    gammacc=0.1128;
    gammaIA=0.1128;
    gammaIc=0.1128;
    gammacAc=0.1128;
    gammaIAc=0.1128;
    deltaA=0.1923;
    deltac=0.1923;
    deltaAc=0.1923;
    eta=1; 
    betaA=0.33296;
    betac=0.23296;
pars = [mu, pi, omegaA,...
    kappaA, kappac, kappaAc,...
    nuA, epsilon, sigmaA,...
    sigmac, sigmaAc, gammacA,...
    gammacc, gammaIA, gammaIc,...
    gammacAc, gammaIAc, deltaA,...
    deltac, deltaAc, eta,...
    betaA, betac];
t0 = 0; 
tf = 50;
% Define the Initial conditions 
S0= 1000;
VA0=400; 
CA0=300;
IA0=100;
RA0=80;
Cc0= 500;
Ic0=200; 
Rc0=100;
CAc0=100; 
IAc0=80; 
RAc0=30;
x0=[S0 VA0 CA0 IA0 RA0 Cc0 Ic0 Rc0 CAc0 IAc0 RAc0];
options = [];
[t,y]=ode45(@Meningitis_coinfection, [t0,tf], x0, options, pars);
figure(1)
plot(t,y(:,1), 'r')
xlabel('time')
ylabel('S')
function dydt = Meningitis_coinfection(t,y,pars)
%Declaring my parameters (pars)
      mu = pars(1); 
      pi = pars(2);   
      omegaA = pars(3);    
      kappaA = pars(4);  
      kappac = pars(5); 
      kappaAc = pars(6); 
      nuA = pars(7);  
      epsilon = pars(8); 
      sigmaA = pars(9);
      sigmac = pars(10);   
      sigmaAc = pars(11); 
      gammacA = pars(12);     
      gammacc = pars(13); 
      gammaIA = pars(14); 
      gammaIc = pars(15);
      gammacAc = pars(16); 
      gammaIAc = pars(17);  
      deltaA = pars(18);
      deltac = pars(19);  
      deltaAc = pars(20); 
      eta = pars(21);  
      betaA = pars(22); 
      betac = pars(23);  
   S =  max(0,y(1));
   VA= max(0,y(2));
   CA = max(0,y(3));
   IA= max(0,y(4));
   RA = max(0,y(5));
   Cc = max(0,y(6));
   Ic= max(0,y(7));
   Rc = max(0,y(8)); 
   CAc= max(0,y(9)); 
   IAc = max(0,y(10));
   RAc = max(0,y(11));
   %Declaring the system of equation (model) 
      N = S + VA + CA + IA + RA + Cc + Ic + Rc + CAc + IAc + RAc ;
      lambdaA = betaA*(eta*CA+IA)/N;   lambdac = betac*(eta*Cc+Ic)/N;
      dS = pi+omegaA*VA+kappaA*RA+kappac*Rc+kappaAc*RAc-(lambdaA+lambdac)*S-(nuA+mu)*S;
      dVA = nuA*S-(1-epsilon)*(lambdaA)*VA-lambdac*VA-(omegaA+mu)*VA;
      dCA = lambdaA*S+(1-epsilon)*lambdaA*VA-lambdac*CA-(sigmaA+gammacA+mu)*CA; 
      dIA = sigmaA*CA-lambdac*IA-(gammaIA+mu+deltaA)*IA; 
      dRA = gammacA*CA+gammaIA*IA-(kappaA+mu)*RA;
      dCc = (S+VA)*lambdac-lambdaA*Cc-(sigmac+gammacc+mu)*Cc; 
      dIc = sigmac*Cc-lambdaA*Ic-(gammaIc+mu+deltac)*Ic; 
      dRc = gammacc*Cc+gammaIc*Ic-(kappac+mu)*Rc;
      dCAc = (Cc+Ic)*lambdaA+lambdac*(CA+IA)-(sigmaAc+gammacAc+mu)*CAc; 
      dIAc = sigmaAc*CAc-(gammaIAc+mu+deltaAc)*IAc; 
      dRAc = gammacAc*CAc+gammaIAc*IAc-(kappaAc+mu)*RAc;
   dydt = [dS;dVA;dCA;dIA;dRA;dCc;dIc;dRc;dCAc;dIAc;dRAc];   
  end
  end
By the way: pi is a constant in Matlab - i recommend not to overwrite.
Best regards
Stephan
Related Question