Help: How to use 'for' loop to plot multiple different values when using while loop for function?

The issue of the code is ' The second 'gr' value not storing in the data to plot the multiple grap on the same figure due to using while loop.

clc; close all; clear all;

%======================SPACEGRID====================================%
ymax=14; m=56; dy=ymax/m; y=dy:dy:ymax;      
xmax=1; n=20;  dx=xmax/n; x=dx:dx:xmax;        
x1=dx:xmax;
%=====================TIMEGRID======================================%
tmax=100; nt=500; dt=tmax/nt; t=dt:dt:tmax;      
%=====================STEADYSTATEINPUTVALUES========================%
tol=1e-2; 
max_difference(1)=1;  max_Iteration=1;
umax_difference(1)=1; %umax_Iteration=1;
%======================TEMPERATUREINPUTVALUES=======================%
pr=6.2; phi=45;
for gr=[10^6,10^8]
%======================EFFECTS=====================================%
M=1; del=-1;
%=======================NFINPUTVALUES===============================%
%for PHI=[0.01 0.4]
PHI=0.01;
rhoF=997.1; rhobetaF=20939.1; rhocpF=4166881; kF=0.613; %H2O
rhoC=8933; rhobetaC=14918.11; rhocpC=3439205; kC=401;   %CU

%==================NF EXPRESSIONS==================================%
KNF=kF*(((kC)+(2*kF)-(2*PHI*(kF-kC)))/((kC)+(2*kF)+(PHI*(kF-kC))));
RHONF=1-PHI+PHI*rhoC/rhoF;
RHOCPNF=1-PHI+PHI*rhocpC/rhocpF;
RHOBETANF=1-PHI+PHI*rhobetaC/rhobetaF;
E1=1/1-PHI^2.5*1/RHONF;
E2=RHOBETANF/RHONF;
E3=1/pr*KNF/kF*1/RHOCPNF;
E4=1/RHONF;
E5=1/RHOCPNF;

%====================INTIALIZATION=================================%
UWALL=0; UOLD=zeros(m,nt); UNEW=0;
VWALL=0; VOLD=zeros(m,nt); VNEW=0;
TNEW=0; TOLD=TNEW*ones(m,nt); TWALL=ones(1,length(t));
A=zeros([1,m]); 
B=A;
C=A;
D=A;
T=TOLD; U=UOLD; V=VOLD;

tic
%===================ENERGYEQUATION==============================%

while max_Iteration>tol
    for j=1:nt
        for i=1:m
            if j>i
                C(i)=dt*VOLD(i,j)/4*dy-dt*E3/(2*dy^2);
            elseif i>j
                A(i)=-dt*VOLD(i,j)/4*dy-dt*E3/(2*dy^2);
            elseif i==j
                B(i)=1+dt*UOLD(i,j)/2*dx+dt*E3/dy^2-dt*E5*del/2;
            end
        end
    end
    for j=2:nt
        for i=2:m-1
            if i==2
                D(i)=-dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/2*dx+dt*E3*(TOLD(i-1,j)-2*TOLD(i,j)+TWALL(j)/(2*dy^2))-dt*VOLD(i,j)*TWALL(j)-TOLD(i-1,j)+TOLD(i,j)/4*dy-dt*(E3+E5)*TWALL(j)/2*dy^2;
            elseif i==m-1
                D(i)=-dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/2*dx+dt*E3*(TNEW-2*TOLD(i,j)+TOLD(i+1,j)/(2*dy^2))-dt*VOLD(i,j)*TOLD(i+1,j)-TNEW+TOLD(i,j)/4*dy-dt*(E3+E5)*TNEW/2*dy^2;
            else
                D(i)=-dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/2*dx+dt*E3*(TOLD(i-1,j)-2*TOLD(i,j)+TOLD(i+1,j)/(2*dy^2))-dt*VOLD(i,j)*TOLD(i+1,j)-TOLD(i-1,j)+TOLD(i,j)/4*dy-dt*E5*del*TOLD(i,j)/2;
            end
        end
        T(:,j)=TriDiag(A,B,C,D);
        dt=0.2+dt;
        T(1,:)=TWALL(j);
        TOLD=T;
%====================STEADY STATE=================================%     
        max_difference(j)=max(abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));
        if max_difference(j)<tol
            break
        end
        max_Iteration=max_difference;
    end
end
%============CALCULATE INTEGRAL BY NUMERICAL===========%
IT=trapz(y,T);

%================CALCULATE PARTIALDERIVATIVE=============%
RT=gradient(IT);%/gradient(x);

%==================MOMENTUM EQUATION=====================%
    for j=1:nt
        for i=1:m
            if j>i
                C(i)=dt*VOLD(i,j)/4*dy-dt*E1/(2*dy^2);
            elseif i>j
                A(i)=-dt*VOLD(i,j)/4*dy-dt*E1/(2*dy^2);
            elseif i==j
                B(i)=1+dt*UOLD(i,j)/2*dx+dt*E1/dy^2+dt*E4*M/2;
            end
        end
    end
    for j=2:nt
        for i=2:m-1
            if i==2
                D(i)=-dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1))/2*dx+dt*E1*(UOLD(i-1,j)-2*UOLD(i,j)+UWALL/2*dy^2)-dt*VOLD(i,j)*(UWALL-UOLD(i-1,j)+UOLD(i,j))/4*dy+dt*E2*cos(phi)*gr^(-1/4)*RT(j)+dt*E2*sin(phi)/2*(TWALL(j)+TOLD(i,j))-dt*(E1+E2+E4)*UWALL/2*dy^2;
            elseif i==m-1
                D(i)=-dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1))/2*dx+dt*E1*(UNEW-2*UOLD(i,j)+UOLD(i+1,j)/2*dy^2)-dt*VOLD(i,j)*(UOLD(i+1,j)-UNEW+UOLD(i,j))/4*dy+dt*E2*cos(phi)*gr^(-1/4)*RT(j)+dt*E2*sin(phi)/2*TNEW+TOLD(i,j)-dt*(E1+E2+E4)*UNEW/2*dy^2;
            else
                D(i)=-dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1))/2*dx+dt*E1*(UOLD(i-1,j)-2*UOLD(i,j)+UOLD(i+1,j)/2*dy^2)-dt*VOLD(i,j)*(UOLD(i+1,j)-UOLD(i-1,j)+UOLD(i,j))/4*dy+dt*E2*cos(phi)*gr^(-1/4)*RT(j)+dt*E2*sin(phi)/2*TNEW+TOLD(i,j)-dt*E4*M*UOLD(i,j)/2;
            end
        end
        U(:,j)=TriDiag(A,B,C,D);
        dt=0.2+dt;
        U(1,:)=UWALL;
        UOLD=U;
%====================STEADYSTATE=================================%          
        umax_difference(j) = max(abs(U(:,j)-U(:,j-1))./max(1,abs(U(:,j-1))));
        if max_difference(j)<tol
            break
        end
        %max_Iteration=max_difference;
    end
%end
figure(1)
if gr==10^6
plot(y,T(:,6),LineWidth=1.5,LineStyle="-",Color='k');
hold on
plot(y,T(:,106),LineWidth=1.5,LineStyle="--",Color='k');
hold on
plot(y,T(:,108),LineWidth=1.5,LineStyle=":",Color='k');
hold on
elseif gr==10^8
plot(y,T(:,6),LineWidth=1.5,LineStyle="-",Color='r');
hold on
plot(y,T(:,106),LineWidth=1.5,LineStyle="--",Color='r');
hold on
plot(y,T(:,108),LineWidth=1.5,LineStyle=":",Color='r');
hold on
end
end
xlim([0.22 14]);ylim([0 4.5]);
legend show
colorbar
toc

Answer 1

[Taylor]

If you're trying to avoid plotting on the same figure, you could move your "figure(1)" command into the "if" statement and place a "figure(2)" into the "elseif" statement. To be completely explicit, you could also call "ax1 = gca;" after "figure(1)" and "ax2 = gca;" after "figure(2)". Is that what you're looking for?