How to plot the graph for different value of t?

I wrote 2D heat equation backward finite difference code. How do I plot the graph for each t?

clc
% Construction of the meshgrid
Lx = 1;
Ly = 1;
Ny = 200;         % Number of grid points in y direction
Nx = 200;         % Number of grid points in x direction
T = 1;
Nt = 3;
dx = Lx/(Nx+1); % Spacing in x direction
dy = Ly/(Ny+1); % Spacing in y direction
dt = T/(Nt+1); % Spacing for the t 
x=0:dx:Lx;
y=0:dy:Ly;
t=0:dt:T;
alpha = 1/(2*pi*pi);
 
[X,Y] = meshgrid(x,y);%2d array of x,y values
X = X' ;             % X(i,j), Y(i,j) are coordinates of (i,j) point
Y = Y';
sx = (alpha*dt)/dx^2;
sy = (alpha*dt)/dy^2;
%-----------------------------------------------------------
Iint = 2:Nx+1; % Indices of interior point in x direction
Jint = 2:Ny+1; % Indices of interior point in y direction
 
Xint = X(Iint,Jint);
Yint = Y(Iint,Jint);
 
%--------Intial condition
u0 = sin(pi*X).*cos(pi*Y);
%uexact = sin(pi*X).*cos(pi*Y).*exp(-1);
 
%------adjust the rhs to include boundary terms
rhs = rhs(u0,Iint, Jint, Nx, Ny ,sx , sy);
 
%-----------------------------------------------
 
rhs;
%convert the rhs into column vector
F = reshape(rhs, Nx*Ny,1);
 
%------assembly of the tridiagonal matrix(LHS)---------
e=ones(Nx,1);
T=spdiags([-sx*e,(2.*sx+2.*sy+1)*e,-sx*e],-1:1,Nx,Nx);
D=spdiags(-sy*e,0,Nx,Nx);
A=blktridiag(T,D,D,Ny);
%-------------------------------------------------
U=zeros(Nx+2,Ny+2,length(t));
for k=1:length(t)
    U(:,:,k)=u0;
end
 
 
for k=2:length(t)
    
    uvec = A\F;
    U(Iint, Jint,k)= reshape(uvec, Nx,Ny);
    rhs = rhs(U,Iint, Jint, Nx, Ny ,sx , sy);
    F = reshape(rhs, Nx*Ny,1);
    
end
 

[dpb]

Which/what style of graph do you want to plot; I've not tried to read the code to know for certain what it returns...

But, I'll note that in the last for...end loop that

 rhs = rhs(U,Iint, Jint, Nx, Ny ,sx , sy);

overwrites the rhs array the first pass through so nothing else will change after that point...

[kaps]

I want to plot 3D plot of the output. I observed this mistake. I replace this as

rhs = rhs(U(:,:,k),Iint, Jint, Nx, Ny ,sx , sy);

function rhs= rhs(u0,Iint, Jint, Nx, Ny ,sx , sy)
rhs = u0(Iint, Jint);
rhs(:,1) = rhs(:,1) + sy.*u0(Iint,1);
rhs(:,Ny)= rhs(:,Ny) + sy.*u0(Iint,Ny+2);
rhs(1,:) = rhs(1,:) + sx.*u0(1,Jint);
rhs(Nx,:) = rhs(Nx,:) + sx.*u0(Nx+2,Jint);
end

[dpb]

That won't fix your problems -- you can't have both a function and a variable of the same name--and I hadn't noticed you had already wiped out the rhs array you started with way early on with the first reassignment.

You need to define that array once and then refer to it elsewhere and not reassign to it anywhere, at least with out subscripting to modify, if needed, certain pieces of it.

Then, if you're going to use a function to do that, name it something different than the variable.

[kaps]

In the actual code, I defined rhs_new function file, then I call it for different U(:,:,k) then I stored this value in rhs variable. Is there any mistake in the following block of code?

rhs = rhs_new(U(:,:,k),Iint, Jint, Nx, Ny ,sx , sy); % calling

%definition 
function rhs= rhs_new(u0,Iint, Jint, Nx, Ny ,sx , sy)
rhs = u0(Iint, Jint);
rhs(:,1) = rhs(:,1) + sy.*u0(Iint,1);
rhs(:,Ny)= rhs(:,Ny) + sy.*u0(Iint,Ny+2);
rhs(1,:) = rhs(1,:) + sx.*u0(1,Jint);
rhs(Nx,:) = rhs(Nx,:) + sx.*u0(Nx+2,Jint);
end