How to make 'for loop' for three variables (i,j,n)?

When I run the code, the result produced same either with the for loop or not. I assumed the for loop is not functioning well in this code. Is there anyone knows how to use 'for loop' for three variables? I want to plot 2D graph

X=zeros(m,p,q);
T=zeros(m,p,q);

%defining boundary conditions
X(:,:,1)=0.35;
X(:,:,q)=0.1;
X(:,1,:)=0.2;

T(:,:,1)=100;
T(:,:,q)=100;
T(:,1,:)=100;

%defining initial condition
X(1,:,:)=0.6;
T(1,:,:)=20;


for i=2:m-1
    for j=2:p-1
        for n=2:q-1

            D=K3*exp(K4*X(i,j,n)-K5/T(i,j,n));
            C=2.3+4.19*X(i,j,n);
            Kr=10.2+79.5*X(i,j,n);
            Kz=5.06+39.8*X(i,j,n);

            part1X=((X(i,j+1,n))^2-2*X(i,j+1,n)*X(i,j-1,n)+(X(i,j-1,n))^2)/(((1-X(i,j,n))^2)*4*delr^2);
            part2X=(X(i,j+1,n)-2*X(i,j,n)+X(i,j-1,n))/((1-X(i,j,n))*delr^2);
            part3X=(X(i,j+1,n)-X(i,j-1,n))/(r(j)*(1-X(i,j,n))*2*delr);
            part4X=(X(i,j,n+1)-2*X(i,j,n)+X(i,j,n-1))/((1-X(i,j,n))*delz^2);
            part5X=((X(i,j,n+1))^2-2*X(i,j,n+1)*X(i,j,n-1)+(X(i,j,n-1))^2)/(((1-X(i,j,n))^2)*4*delz^2);

            part1T=Kr*(T(i,j+1,n)-2*T(i,j,n)+T(i,j-1,n))/delr^2;
            part2T=(T(i,j+1,n)-T(i,j-1,n))*Kr/(r(j)*2*delr);
            part3T=(T(i,j,n+1)-2*T(i,j,n)+T(i,j,n-1))*Kz/delz^2;
            part4T=(X(i+1,j,n)-X(i-1,j,n))*Ps*(C*T(i,j,n)+lamda)/(2*delt);

            X(i+1,j,n)=X(i,j,n)+2*delt*D*(part1X+part2X+part3X+part4X+part5X);
            T(i+1,j,n)=T(i,j,n)+2*delt/((Ps*(0.5+X(i,j,n)+Y)*C)*(part1T-part2T+part3T-part4T));
      
        end   
    end
end

plot(r,X(1,:,2),r,X(2,:,2))
legend("show")
axis([0 2.5 0 0.7])

Answer 1

[Torsten]

for i=1:m-1

instead of

for i=2:m-1

[Ilyani Izzati]

When i=1, it is initial condition (t=0)

[Torsten]

You update the initial data starting with your initial condition:

X(i+1,j,n)=X(i,j,n)+2*delt*D*(part1X+part2X+part3X+part4X+part5X);
T(i+1,j,n)=T(i,j,n)+2*delt/((Ps*(0.5+X(i,j,n)+Y)*C)*(part1T-part2T+part3T-part4T));

Thus the outer loop must start with i = 1, not i = 2.

And all functions have to be evaluated at times with index <= i - otherwise your method had to be implicit.

Thus you will have to correct the setting

part4T=(X(i+1,j,n)-X(i-1,j,n))*Ps*(C*T(i,j,n)+lamda)/(2*delt);

somehow since you work with X(i+1,j,n) that is not yet computed.