Streamline function keeps overwriting plot in a for loop

I am trying to create a streamline plot of a 2D field that is evolving in time. I would like to visualize how the streamlines change over time and I am currently using a for loop to loop through each iteration in time. However, when I use the streamline function to plot it, I can't seem to clear the previous iteration of streamlines and instead they keep layering on top of the previous one. Using the close command on the figure does work, but this causes the figure to flash on the screen for a quick second before it closes again, and I am trying to get it to stay on the screen and keep updating like a normal plot would.

[VBBV]

can you share your code ?

[Shaan Stephen]

I am trying to do the lid driven cavity problem for some CFD code so it is pretty bulky. The section where I print the streamline function is in the fractial step portion

clc; close all; clear;

global np nu nx ny dx dy dt iu iv ip noo uBC_L uBC_R uBC_B uBC_T vBC_L vBC_R vBC_T vBC_B

nx = 64 ;
ny = 64 ;
Re = 1000 ;
noo = 1/Re;

%BCs ...
uBC_L = zeros(1, ny); 
uBC_R = zeros(1, ny);
uBC_B = 0; 
uBC_T = 1;
vBC_L = 0; 
vBC_R = 0; 
vBC_B = 0; 
vBC_T = 0;

% % w/ jet
% midpoint = ny/2;
% midjet = .25;
% uBC_L(midpoint) = midjet; uBC_L(midpoint+1) = midjet; uBC_L(midpoint-1) = midjet;

% Restart?
Restart  = false ;

% pointer and grid
[ip, iu, iv] = GenPointer(nx, ny) ;

% time step

dt   = .01;


dx = 1 / nx ;
dy = 1 / ny ;
np = nx * ny ;
nu = 2*nx*ny - nx - ny ;
nt = 60/dt;

uinitial = zeros(nu, 1);
u_BC_div = BC_Div(uBC_L, uBC_R, vBC_T, vBC_B);
u_BC_Lap = BC_Laplace();


pressure = cell([nt 1]);
% explicit euler
uf = uinitial + dt * (Adv(uinitial) + noo * Laplace(uinitial));
b_EE = (1/dt) * Div(uf) + (1/dt) * u_BC_div;
pressure{1} = CG_EE(b_EE);
u_new = uf - dt * Grad(cell2mat(pressure{1}));

u{1} = uinitial;
u{2} = u_new;
x = 0:dx:1;
y = 0:dy:1;
x = x(2:end-1);
y = y(2:end);


%%
P_array = nan(np, 1);
for i = 1:np
    P_array(i) = i;
end

test1 = Grad(P_array);


h = 0;
xy = cell(nt, 1);
%% time stepping using fractial step
for i = 2:1:nt
    % fractial step: stage 1
    RHS_b = Rinverse(u{i}) + dt * (3 * Adv(u{i})- Adv(u{i-1})) * .5 + dt * noo * u_BC_Lap;
    uf    = CG_stage1(RHS_b) ;
    
    % fractial step: stage 2
    RHS_b = Div(cell2mat(uf)) /dt + BC_Div(uBC_L, uBC_R, vBC_T, vBC_B)/dt;
    pressure{i}  = CG_stage2(RHS_b) ;

    % fractial step: stage 3 (assemble u_new)
    u{i+1} = uf{1} - dt * Grad(cell2mat(pressure{i}));
    s = [nx ny];
    uval = u{end}(1:((nx-1) * ny)); % take only horizontal velocity values
    ugraph = reshape(uval, [ny (nx-1)]);
    %ugraph = [uBC_L ; ugraph; uBC_R];
    vval = u{end}((nx-1)*ny + 1:end);
    vgraph = reshape(vval, [(nx-1) (ny)]);

    % plot
    figure(1)
    contourf(x, y, ugraph)

    % N = 1000;  xstart = max(x)*rand(N,1);  ystart = max(y)*rand(N,1);
    % xstart = max(x)*rand(N,1);  ystart = max(y)*rand(N,1);
    
    % figure(i)
    % xy{i} = stream2(x, y(1:end-1), ugraph(1:end-1, :), vgraph(:, 1:end-1), xstart, ystart,[0.3, 50]);
    % streamlinegraph = streamline(xy{i});
    % 
    i
end

%save('latest_solution.mat', 'u', 'p', ???)

%% testing online code
% close all
uval = u{end}(1:((nx-1) * ny));
vval = u{end}((nx-1)*ny + 1:end);
N = 1000;  xstart = max(x)*rand(N,1);  ystart = max(y)*rand(N,1);
xstart = max(x)*rand(N,1);  ystart = max(y)*rand(N,1);
ugraph = reshape(uval, [ny (nx-1)]);
vgraph = reshape(vval, [(nx-1) ny]);
figure(2)
h = streamline(x, y(1:end-1), ugraph(1:end-1, :), vgraph(:, 1:end-1), xstart, ystart,[0.3, 50])
%% Req'ed functions
% ========================================================================
function [ip, iu, iv] = GenPointer(nx, ny)

% Memory allocation
ip = nan(nx, ny) ;
iu = nan(nx, ny) ;
iv = nan(nx, ny) ;

% Pointer matrix for P
id_p = 1 ; % index to be used in vector variable P
for i = 1:1:nx
    for j = 1:1:ny
        ip(i, j) = id_p ;
        id_p = id_p + 1 ;
    end
end

% Pointer matrix for ux
id_u = 1 ; % index to be used in vector variable u = [ux; uy]
for i = 2:1:nx
    for j = 1:1:ny
        iu(i, j) = id_u ;
        id_u = id_u + 1 ;
    end
end

% Pointer matrix for uy
for i = 1:1:nx
    for j = 2:1:ny
        iv(i, j) = id_u ;
        id_u = id_u + 1 ;
    end
end

end


% ========================================================================
function qo = Grad(qi)
global np nu nx ny dx dy iu iv ip
% Gradient operator: Pressure gradient at u and v points
%       input: p-type (np elements)
%       output: u-type (nu elements)

% Input size check:
if size(qi, 2) ~= 1
    error('Grad:: input is not a column vector')
end

if size(qi, 1) ~= np
    error('Grad:: input vector size mismatch!!!')
end

% Initialize output
qo = nan(nu, 1) ;

% inner domain
% x-direction gradient
for i = 2:1:nx
    for j = 1:1:ny
        qo(iu(i, j)) = ( -qi(ip(i-1, j)) + qi(ip(i, j)) ) / dx ;
    end
end

% y-direction gradient
for i = 1:1:nx
    for j = 2:1:ny
        qo(iv(i, j)) = ( -qi(ip(i, j-1)) + qi(ip(i, j)) ) / dy ;
    end
end

end


% ========================================================================
function qo = Div(qi)
global np nu nx ny dx dy iu iv ip
% divergence operator: To find velocity at p values
%       input: u-type (nu elements)
%       output: p-type (np elements)

% Input size check:
if size(qi, 2) ~= 1
    error('Div:: input is not a column vector')
end

if size(qi, 1) ~= nu
    error('Div:: input vector size mismatch!!!')
end

% Initialize output
qo = nan(np, 1) ;

% inner domain
for i = 2:1:(nx-1)
    for j = 2:1:(ny-1)
        qo(ip(i, j)) = ( - qi(iu(i, j)) + qi(iu(i+1, j)) ) / dx ...
                     + ( - qi(iv(i, j)) + qi(iv(i, j+1)) ) / dy ;
     end
end


% Edges
% bottom inner 
j = 1 ;
for i = 2:1:(nx-1)
    qo(ip(i, j)) = ( - qi(iu(i, j)) + qi(iu(i+1, j)) ) / dx ...
                 + (                  qi(iv(i, j+1)) ) / dy ;   % - qi(iv(i, j)) + 
end
% top inner
j = ny ;
for i = 2:1:(nx-1)
    qo(ip(i, j)) = ( - qi(iu(i, j)) + qi(iu(i+1, j)) ) / dx ...
                 + ( - qi(iv(i, j))                  ) / dy ;   % + qi(iv(i, j+1))
end

% left inner
i = 1 ;
for j = 2:1:(ny-1)
    qo(ip(i, j)) = (                + qi(iu(i+1, j)) ) / dx ... % - qi(iu(i, j)) 
                 + ( - qi(iv(i, j)) + qi(iv(i, j+1)) ) / dy ;
end

% right inner
i = nx ;
for j = 2:1:(ny-1)
    qo(ip(i, j)) = ( - qi(iu(i, j))                  ) / dx ... % + qi(iu(i+1, j)) 
                 + ( - qi(iv(i, j)) + qi(iv(i, j+1)) ) / dy ;
end
% Corners
% bottom left (pinning)
i = 1 ;
j = 1 ;
        qo(ip(i, j)) = 0 ;

% bottom right
i = nx ;
j = 1 ;
        qo(ip(i, j)) = ( - qi(iu(i, j))                  ) / dx ... %  + qi(iu(i+1, j)) 
                     + (                + qi(iv(i, j+1)) ) / dy ;   %  - qi(iv(i, j)) 

% top left
i = 1 ;
j = ny ;
        qo(ip(i, j)) = (                + qi(iu(i+1, j)) ) / dx ... % - qi(iu(i, j)) 
                     + ( - qi(iv(i, j))                  ) / dy ;   % + qi(iv(i, j+1))

% top right
i = nx ;
j = ny ;
        qo(ip(i, j)) = ( - qi(iu(i, j))                  ) / dx ... % + qi(iu(i+1, j)) 
                     + ( - qi(iv(i, j))                  ) / dy ;   % + qi(iv(i, j+1))

end


% ========================================================================
function bcD = BC_Div(uBC_L, uBC_R, vBC_T, vBC_B)
global np ip nx ny dx dy
% BC vector for divergence operator: 
%       input: BCs
%       output: p-type (np elements)

% Initialize output
bcD = zeros(np, 1) ;

% Edges
% bottom inner
j = 1 ;
for i = 2:1:(nx-1)
        bcD(ip(i, j)) = - vBC_B / dy ;   % qi(iv(i, j+1))
end
% top inner
j = ny ;
for i = 2:1:(nx-1)
        bcD(ip(i, j)) =   vBC_T / dy ;   % qi(iv(i, j+1))
end

% left inner
i = 1 ;
for j = 2:1:(ny-1)
    bcD(ip(i, j)) = - uBC_L(j)/ dx ;
end
    
% right inner
i = nx ;
for j = 2:1:(ny-1)
    bcD(ip(i, j)) =   uBC_R(j)/ dx ;
end

% Corners
% bottom left (need pinning)
i = 1 ;
j = 1 ;
    bcD(ip(i, j)) = 0 ;
                 
% bottom right
i = nx ;
j = 1 ;
    bcD(ip(i, j)) = + uBC_R(j) / dx ... 
                    - vBC_B / dy ;         
% top left
i = 1 ;
j = ny ;
    bcD(ip(i, j)) = - uBC_L(j) / dx ... % - qi(iu(i, j)) 
                    + vBC_T / dy ;   % + qi(iv(i, j+1))
             
% top right
i = nx ;
j = ny ;
    bcD(ip(i, j)) = + uBC_R(j) / dx ... 
                    + vBC_T / dy ;   
             
end


% ========================================================================
function qo = Laplace(qi)
global nu iu iv nx ny dx dy 
% Laplace operator: 
%       input: u-type (nu elements)
%       output: u-type (nu elements)

% Input size check:
if size(qi, 2) ~= 1
    error('Div:: Need a colume vector input!!')
end
    
if size(qi, 1) ~= nu
    error('Div:: input vector size mismatch!!!!!')
end

    
% Initialize output
qo = nan(nu, 1) ;

%% 1. ex-Component
% inner domain
for i = 3:1:(nx-1) 
    for j = 2:1:(ny-1)
        qo(iu(i, j)) =   ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ...
                       + ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;
    end
end

% Edges
% left inner 
i = 2 ; 
for j = 2:1:(ny-1)
    qo(iu(i, j)) =   (                 -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
                   + ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;
end

% bottom inner
j = 1 ;
for i = 3:1:(nx-1) 
    qo(iu(i, j)) =   ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ...
                   + ( -qi(iu(i, j))   -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;   % + 2*uBC_B / (dy^2) 
end
      
% right inner
i = nx ; 
for j = 2:1:(ny-1)
    qo(iu(i, j)) =   ( +qi(iu(i-1, j)) -2*qi(iu(i, j))                  ) / (dx^2) ... % + uBC_R / (dx^2)
                   + ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;
end
% top inner
j = ny ;
for i = 3:1:(nx-1) 
    qo(iu(i, j)) =   ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ...
                   + ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) - qi(iu(i, j))   ) / (dy^2) ;   % + 2*uBC_T / (dy^2) 
end
% Corners
% bottom left
i = 2 ; 
j = 1 ;
        qo(iu(i, j)) =   (                 -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ... % + uBC_L   / (dx^2)
                       + ( -qi(iu(i, j))   -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;   % + 2*uBC_B / (dy^2) 

% bottom right
i = nx ; 
j = 1 ;
        qo(iu(i, j)) =   ( +qi(iu(i-1, j)) -2*qi(iu(i, j))                  ) / (dx^2) ... % + uBC_R   / (dx^2)
                       + ( -qi(iu(i, j))   -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;   % + 2*uBC_B / (dy^2) 

% top left
i = 2 ; 
j = ny ;
        qo(iu(i, j)) =   (                 -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ... % + uBC_L   / (dx^2)
                       + ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) - qi(iu(i, j))   ) / (dy^2) ;   % + 2*uBC_T / (dy^2) 

% top right
i = nx ; 
j = ny ;
        qo(iu(i, j)) =   ( +qi(iu(i-1, j)) -2*qi(iu(i, j))                  ) / (dx^2) ... % + uBC_R   / (dx^2)
                       + ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) - qi(iu(i, j))   ) / (dy^2) ;   % + 2*uBC_T / (dy^2) 


%% 2. ey-Component
% inner domain
for i = 2:1:(nx-1) 
    for j = 3:1:(ny-1)
        qo(iv(i, j)) =   ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ...
                       + ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;
    end
end
% Edges
% left inner
i = 1 ; 
for j = 3:1:(ny-1)
    qo(iv(i, j)) =   ( -qi(iv(i, j))   -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + 2*vBC_L / (dx^2)
                   + ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;
end
   
% bottom inner
j = 2 ; 
for i = 2:1:(nx-1)
    qo(iv(i, j)) =   ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... 
                   + (                 -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;   % + vBC_B
end
    
% right inner
i = nx ; 
for j = 3:1:(ny-1)
    qo(iv(i, j)) =   ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) - qi(iv(i, j))   ) / (dx^2) ... % + vBC_R
                   + ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;   % + vBC_B
end
    
% top inner
j = ny ; 
for i = 2:1:(nx-1)
    qo(iv(i, j)) =   ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + vBC_R
                   + ( +qi(iv(i, j-1)) -2*qi(iv(i, j))                  ) / (dy^2) ;   % + vBC_B
end
% Corners
% bottom left
i = 1 ; 
j = 2 ;
        qo(iv(i, j)) =   ( -qi(iv(i, j))   -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + uBC_L   / (dx^2)
                       + (                 -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;   % + 2*uBC_B / (dy^2) 

% bottom right
i = nx ; 
j = 2 ;
        qo(iv(i, j)) =   ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) - qi(iv(i, j))   ) / (dx^2) ... % + uBC_L   / (dx^2)
                       + (                 -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;   % + 2*uBC_B / (dy^2) 

% top left
i = 1 ; 
j = ny ;
        qo(iv(i, j)) =   ( -qi(iv(i, j))   -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + uBC_L   / (dx^2)
                       + ( +qi(iv(i, j-1)) -2*qi(iv(i, j))                  ) / (dy^2) ;   % + 2*uBC_B / (dy^2) 

% top right
i = nx ; 
j = ny ;
        qo(iv(i, j)) =   ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) - qi(iv(i, j))   ) / (dx^2) ... % + uBC_L   / (dx^2)
                       + ( +qi(iv(i, j-1)) -2*qi(iv(i, j))                  ) / (dy^2) ;   % + 2*uBC_B / (dy^2) 

end


% ========================================================================
function bcL = BC_Laplace()
global nu iu iv nx ny dx dy uBC_L uBC_R uBC_B uBC_T vBC_L vBC_R vBC_T vBC_B
% BC vector for divergence operator: 
%       input: BCs
%       output: u-type (nu elements)

% Input size check:
% input BC's are all scalars


% Initialize output
bcL = zeros(nu, 1) ;

%% 1. U-Component
% inner domain


% Edges
% left inner 
i = 2 ; 
for j = 2:1:(ny-1)
    bcL(iu(i, j)) = +   uBC_L(j) / (dx^2) ;
end

% bottom inner
j = 1 ;
for i = 3:1:(nx-1) 
    bcL(iu(i, j)) = + 2*uBC_B / (dy^2) ; 
end

% right inner
i = nx ; 
for j = 2:1:(ny-1)
    bcL(iu(i, j)) = +   uBC_R(j) / (dx^2) ;
end

    
% top inner
j = ny ;
for i = 3:1:(nx-1) 
    bcL(iu(i, j)) = + 2*uBC_T / (dy^2) ; 
end
% Corners
% bottom left
i = 2 ; 
j = 1 ;
        bcL(iu(i, j)) = + uBC_L(j) / (dx^2) + 2*uBC_B / (dy^2) ; 

% bottom right
i = nx ; 
j = 1 ;
        bcL(iu(i, j)) = + uBC_R(j) / (dx^2) + 2*uBC_B / (dy^2) ; 

% top left
i = 2 ; 
j = ny ;
        bcL(iu(i, j)) = + uBC_L(j) / (dx^2) + 2*uBC_T / (dy^2) ; 

% top right
i = nx ; 
j = ny ;
        bcL(iu(i, j)) = + uBC_R(j) / (dx^2) + 2*uBC_T / (dy^2) ; 


%% 2. V-Component
% inner domain

% Edges
% bottom inner 
j = 2 ; 
for i = 2:1:(nx-1)
    bcL(iv(i, j)) = +   vBC_B / (dy^2) ;
end

% top inner
j = ny ; 
for i = 2:1:(nx-1)
    bcL(iv(i, j)) = +   vBC_T / (dy^2) ;
end    

% left inner
i = 1 ; 
for j = 3:1:(ny-1)
    bcL(iv(i, j)) = + 2*vBC_L / (dx^2) ;
end

% right inner
i = nx ; 
for j = 3:1:(ny-1)
    bcL(iv(i, j)) = + 2*vBC_R / (dx^2) ;
end

% Corners
% bottom left
i = 1 ; 
j = 2 ;
        bcL(iv(i, j)) = + 2*vBC_L / (dx^2) + vBC_B / (dy^2) ; 

% bottom right
i = nx ; 
j = 2 ;
        bcL(iv(i, j)) = + 2*vBC_R / (dx^2) + vBC_B / (dy^2) ; 

% top left
i = 1 ; 
j = ny ;
        bcL(iv(i, j)) = + 2*vBC_L / (dx^2) + vBC_T / (dy^2) ; 

% top right
i = nx ; 
j = ny ;
        bcL(iv(i, j)) = + 2*vBC_R / (dx^2) + vBC_T / (dy^2) ; 

end


% ========================================================================
function qo = Adv(qi)

global nu iu iv nx ny dx dy uBC_L uBC_R uBC_B uBC_T vBC_L vBC_R vBC_T vBC_B
% advection operator (BC embedded): -\nabla \cdot (uu) 
%       input: u-type (nu elements)
%       output: u-type (nu elements)
% 

% Input size check:
if size(qi, 2) ~= 1
    error('Adv:: Need a colume vector input!!')
end
    
if size(qi, 1) ~= nu
    error('Adv:: input vector size mismatch!!')
end

% Initialize output
qo = nan(nu, 1) ;

% 1. U-Component
% inner domain
for i = 3:1:(nx-1) 
    for j = 2:1:(ny-1)
        qo(iu(i, j)) = - (1/dx) * ( - ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2 * ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2          ...
                                    + ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 * ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 )   ...
                       - (1/dy) * ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                                    + ( qi(iu(i  ,j  )) + qi(iu(i  ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i  ,j+1)) ) / 2 ) ;             
    end
end    
    
% Edges
% left inner 
i = 2 ;
    for j = 2:1:(ny-1)
        qo(iu(i, j)) = - ( - ( uBC_L(j)        + qi(iu(i  ,j  )) ) / 2 * ( uBC_L(j)        + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 * ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 ) / dx   ...
                       - ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i  ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i  ,j+1)) ) / 2 ) / dy ;                
    end

% bottom inner
j = 1 ;
for i = 3:1:(nx-1)
        qo(iu(i, j)) = - ( - ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2 * ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 * ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 ) / dx   ...
                       - ( -   uBC_B                                   *   vBC_B                                            ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i  ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i  ,j+1)) ) / 2 ) / dy ;         
end 

% right inner
i = nx ;
    for j = 2:1:(ny-1)
        qo(iu(i, j)) = - ( - ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2 * ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + uBC_R(j)        ) / 2 * ( qi(iu(i  ,j  )) + uBC_R(j)        ) / 2 ) / dx   ...
                       - ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i  ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i  ,j+1)) ) / 2 ) / dy ;                
    end

% top inner
j = ny ;
for i = 3:1:(nx-1)
        qo(iu(i, j)) = - ( - ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2 * ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 * ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 ) / dx   ...
                       - ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                           + uBC_T                                     * vBC_T                                     ) / dy ;         
end 

% Corners
% bottom left
i = 2 ;
j = 1 ;
        qo(iu(i, j)) = - ( - ( uBC_L(j)        + qi(iu(i  ,j  )) ) / 2 * ( uBC_L(j)        + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 * ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 ) / dx   ...
                       - ( -   uBC_B                                   *   vBC_B                                            ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i  ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i  ,j+1)) ) / 2 ) / dy ;                


% bottom right
i = nx ;
j = 1 ;
        qo(iu(i, j)) = - ( - ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2 * ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + uBC_R(j)        ) / 2 * ( qi(iu(i  ,j  )) + uBC_R(j)        ) / 2 ) / dx   ...
                       - ( -   uBC_B                                   *   vBC_B                                            ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i  ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i  ,j+1)) ) / 2 ) / dy ;                


% top left
i = 2 ;
j = ny ;
        qo(iu(i, j)) = - ( - ( uBC_L(j)        + qi(iu(i  ,j  )) ) / 2 * ( uBC_L(j)        + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 * ( qi(iu(i  ,j  )) + qi(iu(i+1,j  )) ) / 2 ) / dx   ...
                       - ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                           + uBC_T                                     * vBC_T                                     ) / dy ;         

% top right
i = nx ;
j = ny ;
        qo(iu(i, j)) = - ( - ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2 * ( qi(iu(i-1,j  )) + qi(iu(i  ,j  )) ) / 2          ...
                           + ( qi(iu(i  ,j  )) + uBC_R(j)        ) / 2 * ( qi(iu(i  ,j  )) + uBC_R(j)        ) / 2 ) / dx   ...
                       - ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                           + uBC_T                                     * vBC_T                                     ) / dy ;         


% 2. V-Component
% inner domain
for i = 2:1:(nx-1) 
    for j = 3:1:(ny-1)
        qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                                    + ( qi(iu(i+1,j-1)) + qi(iu(i+1,j  )) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i+1,j  )) ) / 2 )   ...
                       - (1/dy) * ( - ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2 * ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2          ...
                                    + ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 ) ;             
    end
end      
 
% Edges
% left inner 
i = 1 ;
    for j = 3:1:(ny-1)
        qo(iv(i, j)) = - (1/dx) * ( - uBC_L(j)                                  * vBC_L                                         ...
                                    + ( qi(iu(i+1,j-1)) + qi(iu(i+1,j  )) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i+1,j  )) ) / 2 )   ...
                       - (1/dy) * ( - ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2 * ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2     ...
                                    + ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 ) ;             
    end
   
% bottom inner
j = 2 ;
for i = 2:1:(nx-1)     
    qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                                + ( qi(iu(i+1,j-1)) + qi(iu(i+1,j  )) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i+1,j  )) ) / 2 )   ...
                   - (1/dy) * ( - ( vBC_B           + qi(iv(i  ,j  )) ) / 2 * ( vBC_B           + qi(iv(i  ,j  )) ) / 2          ...
                                + ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 ) ;             
end     
% right inner
i = nx ;
    for j = 3:1:(ny-1)
        qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2    ...
                                    + uBC_R(j)                                  * vBC_R                                     )  ...
                       - (1/dy) * ( - ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2 * ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2    ...
                                    + ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 ) ;             
    end
   
% top inner
j = ny ;
for i = 2:1:(nx-1)     
    qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                                + ( qi(iu(i+1,j-1)) + qi(iu(i+1,j  )) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i+1,j  )) ) / 2 )   ...
                   - (1/dy) * ( - ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2 * ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2          ...
                                + ( qi(iv(i  ,j  )) + vBC_T           ) / 2 * ( qi(iv(i  ,j  )) + vBC_T           ) / 2 ) ;             
end  
% Corners
% bottom left
i = 1 ;
j = 2 ;
        qo(iv(i, j)) = - (1/dx) * ( - uBC_L(j)                                  * vBC_L                                         ...
                                    + ( qi(iu(i+1,j-1)) + qi(iu(i+1,j  )) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i+1,j  )) ) / 2 )   ...
                       - (1/dy) * ( - ( vBC_B           + qi(iv(i  ,j  )) ) / 2 * ( vBC_B           + qi(iv(i  ,j  )) ) / 2     ...
                                    + ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 ) ;             

% bottom right
i = nx ;
j = 2  ;
        qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2  ...
                                    + uBC_R(j)                                  * vBC_R                                     )  ...
                       - (1/dy) * ( - ( vBC_B           + qi(iv(i  ,j  )) ) / 2 * ( vBC_B           + qi(iv(i  ,j  )) ) / 2     ...
                                    + ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i  ,j+1)) ) / 2 ) ;             

% top left
j = ny ;
i = 1  ;     
        qo(iv(i, j)) = - (1/dx) * ( - uBC_L(j)                                  * vBC_L                                         ...
                                + ( qi(iu(i+1,j-1)) + qi(iu(i+1,j  )) ) / 2 * ( qi(iv(i  ,j  )) + qi(iv(i+1,j  )) ) / 2 )   ...
                   - (1/dy) * ( - ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2 * ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2          ...
                                + ( qi(iv(i  ,j  )) + vBC_T           ) / 2 * ( qi(iv(i  ,j  )) + vBC_T           ) / 2 ) ;             

% top right
j = ny ;
i = nx ;    
    qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i  ,j-1)) + qi(iu(i  ,j  )) ) / 2 * ( qi(iv(i-1,j  )) + qi(iv(i  ,j  )) ) / 2          ...
                                    + uBC_R(j)                              * vBC_R                                     )  ...
                   - (1/dy) * ( - ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2 * ( qi(iv(i  ,j-1)) + qi(iv(i  ,j  )) ) / 2          ...
                                + ( qi(iv(i  ,j  )) + vBC_T           ) / 2 * ( qi(iv(i  ,j  )) + vBC_T           ) / 2 ) ;             

end
% ========================================================================

function qout = R(qin)
global dt noo
qout = qin - .5 * dt * noo * Laplace(qin);


end
% ========================================================================

function qout = Rinverse(qin)
global dt noo
qout = qin + .5 * dt * noo * Laplace(qin);

end
% ========================================================================


function x_sol = CG_EE(b)
global np
x(:, 1) = rand(np, 1) ;


n = 1;
r(:, n)  = b - Div(Grad(x(:, 1)));
d = r(:, n);
error = 1;
x_sol{n} = x(:, n);
while error > .001
    
    alpha = (r(:, n)' * r(:, n)) / (d' * Div(Grad(d)));
    x(:, n+1) = x(:, n) + alpha * d;
    r(:, n+1) = r(:, n) - alpha * Div(Grad(d));
    beta = (r(:, n+1)' * r(:, n+1))/(r(:, n)' * r(:, n));
    d = r(:, n+1) + beta * d;
    error = norm(r(:, n+1) - r(:, n));
    n = n+1;
    x_sol{n} = x(:, n);

    
end

x_sol = x_sol(end);
end

% ========================================================================

function x_sol = CG_stage1(b)
global nu
x(:, 1) = zeros(nu, 1) ;


n = 1;
r(:, n)  = b - R(x(:, n));
d = r(:, n);
error = 1;
x_sol{n} = x(:, n);
while error > .001
    
    alpha = (r(:, n)' * r(:, n)) / (d' * R(d));
    x(:, n+1) = x(:, n) + alpha * d;
    r(:, n+1) = r(:, n) - alpha * R(d);
    beta = (r(:, n+1)' * r(:, n+1))/(r(:, n)' * r(:, n));
    d = r(:, n+1) + beta * d;
    error = norm(r(:, n+1) - r(:, n));
    n = n+1;
    x_sol{n} = x(:, n);

    
end
x_sol = x_sol(end);
end
% ========================================================================

function x_sol = CG_stage2(b)
global nx ny
x(:, 1) = zeros(nx * ny, 1) ;


n = 1;
r(:, n)  = b - Div(Grad(x(:, n)));
d = r(:, n);
error = 1;
x_sol{n} = x(:, n);
while error > .001
    
    alpha = (r(:, n)' * r(:, n)) / (d' * Div(Grad(d)));
    x(:, n+1) = x(:, n) + alpha * d;
    r(:, n+1) = r(:, n) - alpha * Div(Grad(d));
    beta = (r(:, n+1)' * r(:, n+1))/(r(:, n)' * r(:, n));
    d = r(:, n+1) + beta * d;
    error = norm(r(:, n+1) - r(:, n));
    n = n+1;
    x_sol{n} = x(:, n);

    
end
x_sol = x_sol(end);
end

Answer 1

[Pratyush Swain]

Hi Shaan,

I understand you are unable to clear the the previous iteration of streamlines as it keeps layering on top of the previous ones. One workaround is to use the "cla" command to clear the axes which will delete all visible figures in each iteration. Please refer to the example implementation below:

% Load Dataset %
load wind

% Meshgrid for starting points %
[startX,startY] = meshgrid(80,20:10:50);

% Cell array to hold the streamline object handles %
lineobjs = cell(1,10);

% Looping over new set of data points in each iteration %
for i=1:10
    xi = x(:,:,i);
    yi = y(:,:,i);
    ui = u(:,:,i);
    vi = v(:,:,i);
    
    % Compute the 2-D streamline vertex data %
    verts = stream2(xi,yi,ui,vi,startX,startY);
    
    % Visualize the 2-D matrix of vector fields by calling streamline function %
    lineobjs{i} = streamline(verts);

    % Add title and labels %
    title(['Streamlines at time step ' num2str(i)]);
    xlabel('X');
    ylabel('Y');
    
    % Add pause and then clear the figure object %
    pause(0.5)
    cla 
end

The "cla" command along with "pause" function will help to visualize the streamlines over each iteration without any overlapping.For more information please refer to

1.https://www.mathworks.com/help/matlab/ref/streamline.html

2.https://www.mathworks.com/help/matlab/ref/cla.html

Hope this helps.