please, how can I fix this error; Index in position 1 exceeds array bounds (must not exceed 10001)

Please, I am trying to plot a backward bifurcation diagram but having this errors, kind someone help me.

Errors

Index in position 1 exceeds array bounds (must not exceed 5001).

Error in BIFURCATIONHELLEN (line 50)
while (Root_array(f,1)==0 && Root_array(f,2)==0 && Root_array(f,3)==0)

Please, blow is the code I used. Thank you.

%plotting the bifurcation parameters showing both
%the disease free equilibrium and the endemic equilibrium point.
% to find A, B, C,D for sue in quadratic equation
% File name:BIFURCATION-HELLEN
% saved in Bifur Codes
% parameter values
R0_value=0.5:0.0001:1.5;
Root_array=zeros(length(R0_value),3);
pi=2041;k=0.45;alpha1=1.6;beta=0.35;c=2;alpha2=2;mu=0.02041;
phi=0.06;alpha3=1.3;gamma=0.05;nu=0.2;
d=0.3;alpha4=1.0;omega=0.21;
%eta=0.8;sigma=0.4;
% eta=1.0; mu= 0.055; kappa=10;
% d=0.1; omega=0.5;
% p=0.1; pi=30000; r1=4; gammaf=0.1; gammas =0.00002;
% epsi=0.002; r2=0.08;
% q=0.9;
% sigma=0.5;
hold on
for i=1:1:length(R0_value)
R0=R0_value(i);
R0=(c*beta*(alpha2*k*(gamma+mu)+alpha1*gamma*(1-k)))./(gamma+mu)*(nu+mu+d) ;

Y=(R0*(gamma+mu)*(nu+mu+d))./(c*beta);
p=(nu+mu+d);
A = alpha4*omega*alpha3*phi*pi*(k*alpha2+(1-k)*alpha1);
B= alpha4*omega*pi*Y+alpha3*mu*phi*k*alpha2*pi+alpha3*phi*(1-k)*alpha1*pi*mu+nu*alpha3*phi*k*alpha2*pi+nu*alpha3*phi*(1-k)*...
    alpha1*pi+p*pi*alpha4*omega*phi*alpha3+p*(1-k)*alpha1*pi*...
alpha4*omega+alpha4*omega*nu*pi*phi*alpha3-alpha4*omega*nu*...
(1-k)*alpha1*pi-beta*c*(alpha4*omega*alpha3*phi*k*alpha2*...
pi+alpha4*omega*alpha4*phi*(1-k)*alpha1*pi);
C = mu*pi*Y+nu*pi*Y+p*pi*alpha4*(gamma+mu)*omega+p*pi*...
mu*phi*alpha3+p*(1-k)*alpha1*pi*mu+alpha4*omega*nu*pi*...
(gamma+mu)-beta*c*(pi*alpha4*omega*Y+alpha3*phi*mu*pi*...
(k*alpha2+(1-k)*alpha1));
D = mu*pi*(p*(gamma+mu)-Y);
poly=[A,B,C,D];
r =roots(poly);
len=length(r);
    for t=1:1:len
        if (imag(r(t))~=0) || (real(r(t))<0)
           Root_array(i,t)=0;
        else
Root_array(i,t)=r(t);
        end
    end
end
f=1;
while (Root_array(f,1)==0 && Root_array(f,2)==0 && Root_array(f,3)==0)
        f=f+1;
end
R0_value_Cr=f;
for j=R0_value_Cr:1:length(R0_value)
    Root_array(j,:) =sort(Root_array(j,:));
end
f1=R0_value_Cr;
while (Root_array(f1,2)~=0)
f1=f1+1;
end
R0_value_Cr2=f1;
Zero_1st=R0_value(1,1:R0_value_Cr2-1);
y_zero=zeros(1,length(Zero_1st));
Unstable =R0_value(1,R0_value_Cr:length(R0_value));
plot(Zero_1st,y_zero,'b','LineWidth',3)
plot(Unstable, Root_array(R0_value_Cr:length(R0_value),2),'b','LineWidth',3)
plot(Unstable,Root_array(R0_value_Cr:length(R0_value),3),'r--','LineWidth',3)
xlabel('R_r','FontSize',11)
ylabel('Force of Infection, \lambda^*_r','FontSize',11)
hold off

Answer 1

[Yasasvi Harish Kumar]

Hi,

This is the reason for your error. The condition is always true and f becomes 1002.

f=1;
while (Root_array(f,1)==0 && Root_array(f,2)==0 && Root_array(f,3)==0)
        f=f+1;
end

One correction for this is,

%plotting the bifurcation parameters showing both
%the disease free equilibrium and the endemic equilibrium point.
% to find A, B, C,D for sue in quadratic equation
% File name:BIFURCATION-HELLEN
% saved in Bifur Codes
% parameter values
R0_value=0.5:0.0001:1.5;
Root_array=zeros(length(R0_value),3);
pi=2041;k=0.45;alpha1=1.6;beta=0.35;c=2;alpha2=2;mu=0.02041;
phi=0.06;alpha3=1.3;gamma=0.05;nu=0.2;
d=0.3;alpha4=1.0;omega=0.21;
%eta=0.8;sigma=0.4;
% eta=1.0; mu= 0.055; kappa=10;
% d=0.1; omega=0.5;
% p=0.1; pi=30000; r1=4; gammaf=0.1; gammas =0.00002;
% epsi=0.002; r2=0.08;
% q=0.9;
% sigma=0.5;
hold on
for i=1:1:length(R0_value)
R0=R0_value(i);
R0=(c*beta*(alpha2*k*(gamma+mu)+alpha1*gamma*(1-k)))./(gamma+mu)*(nu+mu+d) ;
Y=(R0*(gamma+mu)*(nu+mu+d))./(c*beta);
p=(nu+mu+d);
A = alpha4*omega*alpha3*phi*pi*(k*alpha2+(1-k)*alpha1);
B= alpha4*omega*pi*Y+alpha3*mu*phi*k*alpha2*pi+alpha3*phi*(1-k)*alpha1*pi*mu+nu*alpha3*phi*k*alpha2*pi+nu*alpha3*phi*(1-k)*...
    alpha1*pi+p*pi*alpha4*omega*phi*alpha3+p*(1-k)*alpha1*pi*...
alpha4*omega+alpha4*omega*nu*pi*phi*alpha3-alpha4*omega*nu*...
(1-k)*alpha1*pi-beta*c*(alpha4*omega*alpha3*phi*k*alpha2*...
pi+alpha4*omega*alpha4*phi*(1-k)*alpha1*pi);
C = mu*pi*Y+nu*pi*Y+p*pi*alpha4*(gamma+mu)*omega+p*pi*...
mu*phi*alpha3+p*(1-k)*alpha1*pi*mu+alpha4*omega*nu*pi*...
(gamma+mu)-beta*c*(pi*alpha4*omega*Y+alpha3*phi*mu*pi*...
(k*alpha2+(1-k)*alpha1));
D = mu*pi*(p*(gamma+mu)-Y);
poly=[A,B,C,D];
r =roots(poly);
len=length(r);
    for t=1:1:len
        if (imag(r(t))~=0) || (real(r(t))<0)
           Root_array(i,t)=0;
        else
Root_array(i,t)=r(t);
        end
    end
end
f=1;
while (Root_array(f,1)==0 && Root_array(f,2)==0 && Root_array(f,3)==0 && f<length(Root_array))
        f=f+1;
end
R0_value_Cr=f;
for j=R0_value_Cr:1:length(R0_value)
    Root_array(j,:) =sort(Root_array(j,:));
end
f1=R0_value_Cr;
while (Root_array(f1,2)~=0)
f1=f1+1;
end
R0_value_Cr2=f1;
Zero_1st=R0_value(1,1:R0_value_Cr2-1);
y_zero=zeros(1,length(Zero_1st));
Unstable =R0_value(1,R0_value_Cr:length(R0_value));
plot(Zero_1st,y_zero,'b','LineWidth',3)
plot(Unstable, Root_array(R0_value_Cr:length(R0_value),2),'b','LineWidth',3)
plot(Unstable,Root_array(R0_value_Cr:length(R0_value),3),'r--','LineWidth',3)
xlabel('R_r','FontSize',11)
ylabel('Force of Infection, \lambda^*_r','FontSize',11)
hold off

Regards