How to make error based on comparison with analytic solution

The end result for this is I am trying to put error norms and the order of the error for each respective norm into a table alongside ndiv. That part isn't the issue. The issue is that the error needs to be a comparison with analytical solution (T_analytic) that updates error for each loop. It also needs to state the number of divisions in the domain required to achieve asymptotic convergence. My code currently gives me NaN for all the error norms and I don't know how to change the code so that error updates by comparison with analytical solution and it outputs error norms correctly at end. Does anyone know how I can do this?

Code

a=25
b=a;
ndiv=[4,8,16,32,64,128,256,512].';
 for nidx=1:length(ndiv)
    nx=ndiv(nidx);
    nz=nx;

x = linspace(0, a, nx);
z = linspace(0, b, nz);
[X, Z] = meshgrid(x,z);
dx = a/(nx-1);
dz = b/(nz-1);


Tnp1 = zeros(nx, nz);
T_analytic=zeros(nx,nz);

Tnp1(:,1) = 20;                                           % Left boundary
Tnp1(:,end) = 20;                                         % Right boundary
Tnp1(1,:) = 20+380*sin((x*pi)/25)+205*sin((x*5*pi)/25);   % Bottom boundary
Tnp1(end,:) = 20;                                         % Top boundary
T_analytic(:,1) = 20;                                           % Left boundary
T_analytic(:,end) = 20;                                         % Right boundary
T_analytic(1,:) = 20+380*sin((x*pi)/25)+205*sin((x*5*pi)/25);   % Bottom boundary
T_analytic(end,:) = 20;   

err = 1;
tol = 1e-8;
k=0;

while  err > tol                                          
   
    Tn = Tnp1;
    k=k+1;
 
    for i = 2:nx-1
        for j = 2:nz-1
            Tnp1(i,j)=Tn(i,j)+1.85*(0.25*(Tn(i+1,j)+Tn(i-1,j)+Tn(i,j+1)+Tn(i,j-1))-Tn(i,j));
           T_analytic(i,j)=20+380/sinh(pi)+205/sinh(5*pi);
        end
    end
    
    err = max(abs(T_analytic(:) - Tnp1(:)));
end
 
All_k(nidx)=k;
   T_numeric = Tnp1;    %Final T Grid for SOR method
% Create a variable for the difference between solutions:
diff = T_numeric(2:end-1,2:end-1) - T_analytic(2:end-1,2:end-1);
% "Unwrap" the 2D array to a 1D row vector;
diffVec = reshape(diff',1,numel(diff));
% Compute Lp norms as usual;
L1norm(nidx,1) = norm(diffVec,1)*(dx*dz);
L2norm(nidx,1) = norm(diffVec,2)*sqrt(dx*dz);
LInfnorm(nidx,1) = norm(diffVec,Inf);



 end
 Norms=table(ndiv,L1norm,L2norm,LInfnorm);

Answer 1

[Walter Roberson]

It also needs to state the number of divisions in the domain required to achieve asymptotic convergence.

Remove the assignment to All_k . Replace it with

Iterations_needed(nidx,1) = k;

And later

Norms = table(ndiv, Iterations_needed, L1norm, L2norm, LInfnorm);

[Walter Roberson]

My code currently gives me NaN for all the error norms

Tnp1(i,j)=Tn(i,j)+1.85*(0.25*(Tn(i+1,j)+Tn(i-1,j)+Tn(i,j+1)+Tn(i,j-1))-Tn(i,j));

For ndiv = 4, those values grow in magnitude without limit, with some of them going to infinity and adjacent ones going to -infinity. Eventually the difference of values of opposite sign becomes an addition that overflows to infinity, and adjacent entries overflow -infinity. Then the sum or difference of infinities gives you Nan. err then becomes nan, and so err > something is false because comparing anything to nan is false.