Double For cycle plotting

Hi, i need help about plotting those double for cycle in this code. When i plot the 2 for cycle, at the interface of the multilayer soil (at Z1 and T1) it seems that matlab doesn't see the 2°layer from Z1 to Z2 and from T1 to T2. The code should represent the trend of the concentration of a solute as a function of depth and time, like the image in the bottom of this question.

those images are old but it should appear something like this. The problem is that i cant connect the plot of the 2 for cycle in the same graph. Please help me <3

function twolayerupwind_2

clc
clear all

%This is the fucking file of the exam exercise
%Explicit Euler Method is used for solving the following
%initial-boundary-value problem:
%Ut+v1*Ux=alpha1Uxx     Ut+v2*Ux=alpha2Uxx
%U(z,0)=0                                                                      %Cond iniziale
%U(0,t)=g1(t)       U(L,t)=g2(t)                                               %Cond. Al bordo

%Problem parameters

%---------------------layer 1----------------------
alpha1 =   20 ; th1=1; Z1 = 100; %T = 10; -
nx = 300; nt1 = 500; v1=40; Dm1=alpha1; C0=50; T1=Z1/v1;
%---------------------Interface--------------------
xh=Z1;
                                                                        %profondità dell'interfaccia.
%---------------------layer 2----------------------
alpha2 =  2 ;  Z2 = 100;  v2=5 ; Dm2=alpha2; th2=th1*(v1/v2); 
Z=Z1+Z2; nt2 = nt1; T2=Z2/v2;
%Check of problem parameters
if any([alpha1 alpha2  Z1 Z2  T1 T2 Z-xh nx-2 nt1-2 nt2-2]<=0)
    error('Check problem parameters')
end



%Stability Layer 1
st = 1; dx = Z/nx;
while 2*alpha1*T1/nt1/dx^2 + T1/nt1/dx*abs(v1) > st
nt1=nt1+1;
end

%Stability Layer 2
while 2*alpha2*T2/nt2/dx^2 + T2/nt2/dx*abs(v2) > st
nt2=nt2+1;
end


%Post-stabilità andiamo allora ad esplicitare i seguenti termini:
dt1 = T1/nt1; r1 = dt1/dx^2; s1 = dt1/dx;
dt2 = T2/nt2; r2 = dt2/dx^2; s2 = dt2/dx;


%Inizialization
nt1= (T1/dt1)+1
nt2= (T2/dt2)+1

h= nx*xh/Z + 1; h = double(uint16(abs(h))); 


x=linspace(0,Z,nx+1); t1=linspace(0,T1,nt1+1); t2=linspace(0,T2,nt2+1);
xh = x(h);
u=getPhi(x,nx,Z); U=u;                                                     %Fissiamo la nostra funzione concentrazione(u) pari prorpio alla

if v1>=0
p1 = alpha1*r1 + s1*v1; q1 = alpha1*r1;
else
p1 = alpha1*r1; q1 = alpha1*r1 - s1*v1;
end

if v2>=0
p2 = alpha2*r2 + s2*v2; q2 = alpha2*r2;
else
p2 = alpha2*r2; q2 = alpha2*r2 - s2*v2;
end


%Subfunction 
%[Per le condizioni iniziali]
    function Pisciazz = getPhi(x,nx,Z)
        Pisciazz(1:nx+1,1)=0;
        %Pisciazz(1:h,1)=0                                                               %Cond iniziale asse z, 1°Strato
        %Pisciazz(h+1:nz+1,1)=0                                                          %Cond iniziale asse z, 2°Strato
    end
% [Per bordo z=0] 
function VarTempo = getVarTempo(t1)
%VarTempo=2
Dt=0.5; C0=50;
if t1<=Dt
    VarTempo= C0;
else
    VarTempo= 0;
end

end

%Subfunction [Per bordo z=L]
function NoemannNothing = getNoemannNothing(t)
NoemannNothing= u(end-1);

end

%Finite Difference Method [Forward for the time and Central for the space]

for j=2:nt1+1
    u(1)=getVarTempo(t1(j));
    u(2:h-1)= p1*u(1:h-2) + (1- p1 - q1 )*u(2:h-1) + q1*u(3:h);
    %buona u(h)=(th1*Dm1*(-4*u(h-1)+u(h-2))-th2*Dm2*((4*u(h+1)-u(h+2))))/(-3*(th1*Dm1+th2*Dm2));
    %u(h)=(th1*Dm1*(-4*u(h-1)+u(h-2))-th2*Dm2*((4*u(h+1)-u(h+2))))/(((v1*th1-v2*th2)*2*dx)-3*(th1*Dm1+th2*Dm2));
    u(h)=u(h-1);
    
    figure(1);
    plot(u,x,'r*:',U,x,'k');
    xlabel('C(z,t1)'); ylabel('z');
    title(['t=',num2str(t1(j)),'seconds']);
    axis([0 50 0 Z1+Z2]);
    pause(0.01);
    A(j) = trapz(x,u)  
    
end

for j=nt1+2:nt2+1
    u(h)=u(h-1);
    u(h+1:end-1)= p2*u(h:end-2) + (1- p2 - q2 )*u(h+1:end-1) + q2*u(h+2:end);
    u(end)=getNoemannNothing(t);

    figure(1);
    plot(u,x,'r*:',U,x,'k');
    xlabel('C(z,t2)'); ylabel('z');
    title(['t=',num2str(t2(j)),'seconds']);
    axis([0 50 0 Z1+Z2]);
    pause(0.01);
    A(j) = trapz(x,u)  
end


legend(['t=',num2str(T1)],'t=0');
figure(2)
plot(0:dt1:T1,A)

end

Answer 1

[Nehemiae]

Hello,

In order to make the two for loops update the same figure, the “hold on” command can be used in MATLAB. The code written below, and the changes made to ensure that each for loop updates the figure each iteration and clears the unnecessary points, while ensuring that both for loops affect the same figure. Also, I have made changes in the iteration conditions of the second for loop, or else it never executes.

figure; 

for j=2:nt1+1 
    u(1)=getVarTempo(t1(j)); 
    u(2:h-1)= p1*u(1:h-2) + (1- p1 - q1 )*u(2:h-1) + q1*u(3:h); 
    u(h)=u(h-1); 
    
    plot(u,x,'r*:',U,x,'k'); 
    xlabel('C(z,t1)'); ylabel('z'); 
    title(['t=',num2str(t1(j)),'seconds']); 
    axis([0 50 0 Z1+Z2]); 
    pause(0.01); 
    A(j) = trapz(x,u);   
end 

hold on 

for j=nt1+1:2*nt2 
    u(h)=u(h-1); 
    u(h+1:end-1)= p2*u(h:end-2) + (1- p2 - q2 )*u(h+1:end-1) + q2*u(h+2:end); 
    u(end)=getNoemannNothing(t1); 
    
    plot(u,x,'r*:',U,x,'k'); 
    xlabel('C(z,t2)'); ylabel('z'); 
    title(['t=',num2str(t2(j-nt1)),'seconds']); 
    axis([0 50 0 Z1+Z2]); 
    pause(0.01); 
    A(j) = trapz(x,u); 
    hold off 
end 

The resulting graph is as follows.

The documentation on the hold command (https://www.mathworks.com/help/matlab/ref/hold.html) is useful in understanding the above code.