No obvious error in your output? I do not see any error message ?
As outside observers who have not been given documentation, we must assume that this is what your code was designed to calculate.
Otherwise we have to start making random guesses like suggesting that your r should be different, or that your f3 function should involve x4 instead of x3... because how would we know?
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/822230/Untitled.m';
S = urlread(filename)
S =
'clc;
close all;
clear all;
S0=100;E0=10;I0=3;R0=0;D0=0;
lambda=0;a1=0.3;a2=0.4;r=0.23;K=0.3;phi=0.25;b=0.45;mu=0.12;
beta=1;h=0.02;N=5;
tfinal=100;
f1=@(x1,x2,x3,x4,x5) lambda-a1*x1*x2-a2*x1*x3+r*x4;
f2=@(x1,x2,x3,x4,x5) a1*x1*x2+a2*x1*x3-K*x2-phi*x2;
f3=@(x1,x2,x3,x4,x5) K*x2-(b+mu)*x3;
f4=@(x1,x2,x3,x4,x5) b*x3+phi*x2-r*x4;
f5=@(x1,x2,x3,x4,x5) mu*x3;
S(N+1)=[0];E(N+1)=[0];I(N+1)=[0];R(N+1)=[0];D(N+1)=[0];
Sp(N+1)=[0];Ep(N+1)=[0];Ip(N+1)=[0];Rp(N+1)=[0];Dp(N+1)=[0];
M1=0;M2=0;M3=0;M4=0;M5=0;L1=0;L2=0;L3=0;L4=0;L5=0;
t=0:tfinal/N:tfinal;
Sp(1)=S0+h^(beta)*(f1(S0,E0,I0,R0,D0))/gamma(beta)*(beta);
Ep(1)=E0+h^(beta)*(f2(S0,E0,I0,R0,D0))/gamma(beta)*(beta);
Ip(1)=I0+h^(beta)*(f3(S0,E0,I0,R0,D0))/gamma(beta)*(beta);
Rp(1)=R0+h^(beta)*(f4(S0,E0,I0,R0,D0))/gamma(beta)*(beta);
Dp(1)=D0+h^(beta)*(f5(S0,E0,I0,R0,D0))/gamma(beta)*(beta);
S(1)=S0+h^(beta)*(f1(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f1(S0,E0,I0,R0,D0))/gamma(beta+2);
E(1)=E0+h^(beta)*(f2(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f2(S0,E0,I0,R0,D0))/gamma(beta+2);
I(1)=I0+h^(beta)*(f3(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f3(S0,E0,I0,R0,D0))/gamma(beta+2);
R(1)=R0+h^(beta)*(f4(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f4(S0,E0,I0,R0,D0))/gamma(beta+2);
D(1)=D0+h^(beta)*(f5(Sp(1),Ep(1),Ip(1),Rp(1),Dp(1)))/gamma(beta+2)+h^(beta)*beta*(f5(S0,E0,I0,R0,D0))/gamma(beta+2);
for n=1:N
M1=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f1(S0,E0,I0,R0,D0));
M2=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f2(S0,E0,I0,R0,D0));
M3=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f3(S0,E0,I0,R0,D0));
M4=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f4(S0,E0,I0,R0,D0));
M5=((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1))*(f5(S0,E0,I0,R0,D0));
L1=((n+1)^(beta)-(n)^(beta))*(f1(S0,E0,I0,R0,D0));
L2=((n+1)^(beta)-(n)^(beta))*(f2(S0,E0,I0,R0,D0));
L3=((n+1)^(beta)-(n)^(beta))*(f3(S0,E0,I0,R0,D0));
L4=((n+1)^(beta)-(n)^(beta))*(f4(S0,E0,I0,R0,D0));
L5=((n+1)^(beta)-(n)^(beta))*(f5(S0,E0,I0,R0,D0));
for j=1:n
M1=M1+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f1(S(j),E(j),I(j),R(j),D(j)));
M2=M2+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f2(S(j),E(j),I(j),R(j),D(j)));
M3=M3+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f3(S(j),E(j),I(j),R(j),D(j)));
M4=M4+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f4(S(j),E(j),I(j),R(j),D(j)));
M5=M5+((n-j+2)^(beta+1)+(n-j)^(beta+1)-2*(n-j+1)^(beta+1))*(f5(S(j),E(j),I(j),R(j),D(j)));
L1=L1+((n+1-j)^(beta)-(n-j)^(beta))*(f1(S(j),E(j),I(j),R(j),D(j)));
L2=L2+((n+1-j)^(beta)-(n-j)^(beta))*(f2(S(j),E(j),I(j),R(j),D(j)));
L3=L3+((n+1-j)^(beta)-(n-j)^(beta))*(f3(S(j),E(j),I(j),R(j),D(j)));
L4=L4+((n+1-j)^(beta)-(n-j)^(beta))*(f4(S(j),E(j),I(j),R(j),D(j)));
L5=L5+((n+1-j)^(beta)-(n-j)^(beta))*(f5(S(j),E(j),I(j),R(j),D(j)));
end
Sp(n+1)=(S0+(1/gamma(beta))*L1);
Ep(n+1)=(E0+(1/gamma(beta))*L2);
Ip(n+1)=(I0+(1/gamma(beta))*L3);
Rp(n+1)=(R0+(1/gamma(beta))*L4);
Dp(n+1)=(D0+(1/gamma(beta))*L5);
S(n+1)=S0+(h)^(beta)*((f1(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M1)/gamma(beta+2);
E(n+1)=E0+(h)^(beta)*((f2(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M2)/gamma(beta+2);
I(n+1)=I0+(h)^(beta)*((f3(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M3)/gamma(beta+2);
R(n+1)=R0+(h)^(beta)*((f4(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M4)/gamma(beta+2);
D(n+1)=D0+(h)^(beta)*((f5(Sp(n+1),Ep(n+1),Ip(n+1),Rp(n+1),Dp(n+1)))+M5)/gamma(beta+2);
end
Y1(1)=S0;
Y2(1)=E0;
Y3(1)=I0;
Y4(1)=R0;
Y5(1)=D0;
for k=1:N
Y1(k+1)=S(k);
Y2(k+1)=E(k);
Y3(k+1)=I(k);
Y4(k+1)=R(k);
Y5(k+1)=D(k);
end
plot(t,Y1)
'
eval(S)