How to plot 2D line graph to compare the approximate solution with the actual solution?

I am trying to solve the following problem.

I wrote the following code to solve this problem. I stored my approximate solution in the matrix U of size Nx+2 times Ny+2. I want to compare my approximated solution with the actual solution. I am not sure how to plot a matrix in to 2D line graph. Can anyone please help me with this step?

%------Construction of meshgrids-------------------
Lx = 1;
Ly = 1;
%Ny=Nx;
Nx = 3;     % Number of interior grid points in x direction
Ny = 3;     % Number of interior grid points in y direction
dx = Lx/(Nx+1); % Spacing in the x direction
dy = Ly/(Ny+1); % Spacing in the y direction
x=0:dx:Lx;
y=0:dy:Ly;

[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';
% -------------------------------------------------------
Iint = 1:Nx+1; % Indices of interior point in x direction
Jint = 1:Ny+2; % Indices of interior point in y direction
Xint = X(Iint,Jint);
Yint = Y(Iint,Jint);
U = zeros(Nx+2,Ny+2); % Define U to store the answer
%--------------------------------------------------------
uinit = zeros(Nx+1,Ny+2);
u0 = @(x,y) cos(2*pi*x).*cos(2*pi*y);
U(1,:) = u0(x(1),y(Jint));
U(Nx+2,:)=U(1,:);
uinit(1,:)=U(1,:);
uinit(Nx,:)=uinit(1,:);
F1 = reshape(uinit, (Nx+1)*(Ny+2),1);
%---------------------------------------------
%assembly of the tridiagonal matrix(LHS)

sx = 1/(dx^2);
sy = 1/(dy^2);

e=ones(Nx,1);
A = zeros(Nx,Nx+2);
B = zeros(Nx+1,Nx+2);
T=spdiags([sx.*e,(-2*sx)+(-2*sy).*e,sx.*e],-1:1,Nx,Nx);
A(:,2:Nx+1) = T;
B(1:Nx,:)=A;
B(Nx+1,1)=-1;
B(Nx+1,2)=1;
B(Nx+1,Nx+1)=1;
B(Nx+1,Nx+2)=-1;


%T(1,Nx+1)= sx;
%T(Nx+1,1)= sx;
D = zeros(Nx+1,Nx+2);
D1 = zeros(Nx,Nx+2);
D2=spdiags(sy.*e,0,Nx,Nx);
D1(:,2:Nx+1)=D2;
D(1:Nx,:)=D1;
C=blktridiag(B,D,D,Ny+2);


for i=1:Nx
    for j=1:Nx+1
        C(i,j)=(1/2).*C(i,j);
        C((Nx+1)*(Ny+1)+i,(Nx+2)*(Ny+1)+j)=(1/2).*C((Nx+1)*(Ny+1)+i,(Nx+2)*(Ny+1)+j);
    end
end

%---------------------------------------------------------------
%----------Compuating the R.H.s term----------------------------
f = @(x,y) (-8*pi*pi).*cos(2*pi*x).*cos(2*pi*y);
%Solve the linear system
rhs = f(Xint,Yint);
rhs_new = zeros(Nx+1,Ny+2);
rhs_new(1:Nx,:)=rhs(2:Nx+1,:);

%for j=1:Ny+2
%    for i=2:Nx+3
%        rhs(i,j) = (-8*pi*pi)*cos(2*pi*x(i)).*cos(2*pi*y(j));
%    end
%end

for i=1:Nx
    rhs_new(i,1)=(1/2).*rhs_new(i,1);
    rhs_new(i,Ny+2)=(1/2).*rhs_new(i,Ny+2);
end

%convert the rhs into column vector
F = reshape(rhs_new, (Nx+1)*(Ny+2),1);
F2 = F - sx.* F1;




uvec = C\F2;
v= reshape(uvec, Nx+2,Ny+2);
U(2:Nx+1,:)=v(2:Nx+1,:);



%----------Computation of Exact solution--------------------
ue = zeros(Nx+2,Ny+2);
ue(:,:) = cos(2*pi*X) .* cos(2*pi*Y); % Exact Solution

%Error
err = norm(U(:,:) - ue(:,:))/norm(ue(:,:));