Dear all

I have simulated the CTMC model while facing an issue of array bound error

anyway who can resolve the error using the code here attached

%CTCM for transmission dynamics of rabies

%By Mfano Charles

%PhD Student (AMCS) - NM-AIST, Arusha, Tz

%21 June, 2024

clear all

close all

clc;

set(0,'DefaultAxesFontSize', 10)%Increases axes labels

set(gca,'fontsize',18);%Changes font size in a figure.

%-------------------------------------------------------------------------

%% Parameter values of the model

time=79; %time

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%

rho_1=10; %Deterrent coeffiecient of domestic dogs from I_F

rho_2=8; %Deterrent coeffiecient of domestic dogs from I_D

rho_3=15; %Deterrent coeffiecient of domestic dogs from M

theta_1=2000; % recurement of susceptible human

tau_1=0.0004; %The rate that S_H get infections from I_F

tau_2=0.0006; %The rate S_H get infection from I_D %paremeter Estimation

tau_3=0.0003; %The rate that S_H get infection from M

beta_1=1/6; %progression rate out of E_H to I_H(Incubation period of human)

beta_2=1; %Recovery rate of E_H

beta_3=0.54; %Rate of immunity loss of human being

mu_1=0.0142; %Natural death rate of human

sigma_1=1; %Disease induced rate of I_H

theta_2=1000; %Recruitment rate of free range dogs

gamma=1/6; %Progression rate of E_F to I_F

kappa_1=0.00002; %The rate that S_F becomes I_F from I_F

kappa_2=0.00008;%The rate that S_F bceomes infection from I_D

kappa_3=0.00004; %The rate that S_F gets infection form M

sigma_2=0.09;% Disease induced death rate of I_F

mu_2=0.067; %Natural Mortality rate of free-range dogs

theta_3=1200;%Recruitment rate of domestic dog population

psi_1=0.00009; %The rate that S_D get infection from I_F

psi_2=0.00007; %The rate that S_D get infection from I_D

psi_3=0.00003; %The rate that S_D get infection from M

sigma_3= 0.09;%1; %Disease induced death rate for I_D

gamma_1=1/6; %The rate at which E_D becomes I_D

gamma_2=0.09; %Recovery rate of E_D

gamma_3=0.05; %The rate of loss of immunity of E_D

mu_3=0.08; % The naural Mortality rate of domestic dogs

nu_1=0.002; %Environment shedding rate from I_H

nu_2=0.009; %Environment shedding rate from I_F

nu_3=0.005; %Environment shedding rate from I_D

mu_4=0.08; %Naural removal rate of rabies

C=0.003; % Rabies concetration in the environment

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%

%% Euler’s method for solving ODE

S_H= theta_1/mu_1;

S_F= theta_2/mu_2;

S_D= theta_3/mu_3;

%++++++++++++++++++++++++++++++++++++%

%X=S_H; X1=E_H; X2=I_H; X3=R_H; P=S_F; P1=E_F; P2=I_F; Y=S_D;Y1=E_D;Y2=I_D;

%Y3=R_D; Z=M

%+++++++++++++++++++++++++++++++++++++++++++++++++%

%x=S_H; x1=E_H; x2=I_H; x3=R_H; p=S_F; p1=E_F; p2=I_F; y=S_D;y1=E_D;y2=I_D;

%y3=R_D; z=M

%++++++++++++++++++++++++++++++++++++%

dt=0.05; tim=time/dt; kk=0:dt:time;

%Initial conditions

a0=1; b0=0;c0=1; d0=1; e0=1;f0=1;g0=1;q0=1; k0=1; h0=2;

m0=1;n0=1;

%---------------------------------------------------------------------------

%x(1)=a0;

x1(1)=b0;

x2(1)=c0;

x3(1)=d0;

%p(1)=e0;

p1(1)=f0;

p2(1)=g0;

%y(1)=q0;

y1(1)=k0;

y2(1)=h0;

y3(1)=m0;

z(1)=n0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x(1)=S_H-(b0+c0+d0);

p(1)=S_F-(f0+g0);

y(1)=S_D-(k0+h0+m0);

%--------------------------------------------------------------------------------

for jj=1:tim

% HUMAN POPULATION

x(jj+1) = x(jj) + dt * (theta_1 + beta_3 * x3(jj) - mu_1 * x(jj) - (tau_1 * p2(jj) + tau_2 * y2(jj) + tau_3 * z(jj) / (z(jj) + C)) * x(jj));

x1(jj+1) = x1(jj) + dt * ((tau_1 * p2(jj) + tau_2 * y2(jj) + tau_3 * z(jj) / (z(jj) + C)) * x(jj) - (mu_1 + beta_1 + beta_2) * x1(jj));

x2(jj+1) = x2(jj) + dt * (beta_1 * x1(jj) - (sigma_1 + mu_1) * x2(jj));

x3(jj+1) = x3(jj) + dt * (beta_2 * x1(jj) - (beta_3 + mu_1) * x3(jj));

% FREE RANGE POPULATION

p(jj+1) = p(jj) + dt * (theta_2 - (kappa_1 * p2(jj) + kappa_2 * y2(jj) + kappa_3 * z(jj) / (z(jj) + C)) * p(jj) - mu_2 * p(jj));

p1(jj+1) = p1(jj) + dt * ((kappa_1 * p2(jj) + kappa_2 * y2(jj) + kappa_3 * z(jj) / (z(jj) + C)) * p(jj) - (mu_2 + gamma) * p1(jj));

p2(jj+1) = p2(jj) + dt * (gamma * p1(jj) - (mu_2 + sigma_2) * p2(jj));

% DOMESTIC DOG POPULATION

y(jj+1) = y(jj) + dt * (theta_3 - ((psi_1 / (1 + rho_1)) * p2(jj) + (psi_2 / (1 + rho_2)) * y2(jj) + (psi_3 * z(jj) / (1 + rho_3)) * (z(jj) + C)) * y(jj) - mu_3 * y(jj) + gamma_3 * y3(jj));

y1(jj+1) = y1(jj) + dt * (((psi_1 / (1 + rho_1)) * p2(jj) + (psi_2 / (1 + rho_2)) * y2(jj) + (psi_3 * z(jj) / (1 + rho_3)) * (z(jj) + C)) * y(jj) - (mu_3 + gamma_1 + gamma_2) * y1(jj));

y2(jj+1) = y2(jj) + dt * (gamma_1 * y1(jj) - (mu_3 + sigma_3) * y2(jj));

y3(jj+1) = y3(jj) + dt * (gamma_2 * y1(jj) - (mu_3 + sigma_3) * y3(jj));

% RABIES IN THE ENVIRONMENT

z(jj+1) = z(jj) + dt * ((nu_1 * x2(jj) + nu_2 * p2(jj) + nu_3 * y2(jj)) - mu_4 * z(jj));

%================================================================

end

%Three sample paths for CTMC (Gillespie’s algorithm)

for k=1:3

clear t X X1 X2 X3 P P1 P2 Y Y1 Y2 Y3 Z

j=1;

% X1(1)=a0; X2(1)=b0; P1(1)=c0; P2(1)=m0;%P2(1)=w*c0;

% Y1(1)=k0; Y2(1)=n0;%Y2(1)=eta*k0;

% Z(1)=h0;

%X(1)=a0;

X1(1)=b0;

X2(1)=c0;

X3(1)=d0;

%P(1)=e0;

P1(1)=f0;

P2(1)=g0;

%Y(1)=q0;

Y1(1)=k0;

Y2(1)=h0;

Y3(1)=m0;

Z(1)=n0;

X(1)=S_H-(b0+c0+d0);Y(1)=S_F-(e0+f0);P=S_D-(k0+h0+m0);

t(1)=0; %Starting at zero

%------------------------------------------------------------

while (X1(j)+X2(j)+P1(j)+P2(j)+Y1(j)+Y2(j)+Z(j))>0 && t(j)