Please help me to find the errors in my code for predictor corrector method.

I have written the code but it is unable to produce correct result. Can anybody help me to find the errors in the above code?

[Walter Roberson]

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)

Answer 1

[Jan]

I do not know, what you consider as "correct" result. All I see is the running code and there is no chance to predict, what you wnt to do instead.

But at least I can simplify your code, such that a debugging is much easier. Instead of repeating any equation 5 times, I've converted [f1,f2, f3, f4, f5] to one function f and did the same for the 5 inputs:

function Untitled2

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;
t      = 0:tfinal/N:tfinal;

f = @(x) [ ...
   lambda - a1 * x(1) * x(2) - a2 * x(1) * x(3) + r * x(4); ...
   a1 * x(1) * x(2) + a2 * x(1) * x(3) - K * x(2) - phi * x(2); ...
   K * x(2) - (b + mu) * x(3); ...
   b * x(3) + phi * x(2) - r * x(4); ...
   mu * x(3)];

S  = zeros(5, N+1);
Sp = zeros(5, N+1);

S0       = [100; 10; 3; 0; 0];
Sp(:, 1) = S0 + h^beta * f(S0) / gamma(beta) * beta;
S(:, 1)  = S0 + h^beta * f(Sp) / gamma(beta+2) + h^beta * beta *...
   f(S0) / gamma(beta+2);
for n=1:N
   M = ((n+2)^(beta+1)+n^(beta+1)-(2*(n+1))^(beta+1)) * f(S0);
   L = ((n+1)^beta - n^beta) * f(S0);
   
   for j=1:n
      M = M + ((n-j+2)^(beta+1) + (n-j)^(beta+1) - 2*(n-j+1)^(beta+1)) * ...
         f(S(:, j));      
      L = L + ((n+1-j)^beta - (n-j)^beta) * f(S(:, j));
   end
   Sp(:, n+1) = S0 + L / gamma(beta);
   S(:, n+1) = S0 + h^(beta) * (f(Sp(:, n+1)) + M) / gamma(beta+2);
end
Y1 = [S0, S(:, 1:N)];
plot(t, Y1)
end

Adding spaces around the operators improve the readability.

It looks strange, that you calculate S(:, 6), but do not use it in the output.

[Surath Ghosh]

Thank you very much for your comment. Actually Sp, Ip, Ep...... are the predicted values and S, I, E....are the corrected values. I am getting wrong result. I do not understand about the errors in my logic. Can you explain me in details? This is the code of a numerical method that is predictor - Corrector method in mathematics.

[Surath Ghosh]

That is okay. Thank you very much for simplifying. Actually, I am varyfying a research paper. According to that paper values of S will be always less than 100 in the interval (0, 100). Here, I have just plotted S(t) for checking. I have attached the graph. Please check it.