Answers keep on iterating

This program is made to approximate the solution to the Poisson equation. I get the correct answer but it keeps on repeating the same answers (I think 12 times). I only need one set of these answers (as highlighted). Please help, thank you so much in advance.

%Initialization

a=0;
b=2;
c=0;
d=1;
m=5;
n=6;
N= 100;
TOL = 10^-10;

syms ga(y) gb(y) gc(x) gd(x) f(x,y);
f(x,y)=x*exp(y);
ga(y)=0;  %g(a,y)
gb(y)=2*exp(y);  %g(b,y)
gc(x)=x;          %g(x,c)           
gd(x) =x*exp(1);   %g(x,d)

%Step 1
h=(b-a)/n;
k=(d-c)/m;

%Step 2;
for i = 1:n-1
    xi(i)=a+i*h;   %x sub i (k) = xi(k)
end

%Step 3
for j = 1:m-1
    yj(j)=c+j*k;
end

%Step 4
for i=1:n-1
    for j=1:m-1
        wij(i,j)=0;
    end
end

%Step 5
lambda=h^2/k^2;
mu=2*(1+lambda);
l=1;

%Step 6
while l <= N
    
    %Step 7
    z=(-(h^2)*(f(xi(1),yj(m-1)))+ga(yj(m-1))+lambda*gd(xi(1))+lambda*wij(1,m-2)+wij(2,m-1))/mu;
    NORM=abs(z-wij(1,m-1));
    wij(1,m-1)=z;
    
    %Step 8
    for i=2:n-2
        z = (-h^2*f(xi(i),yj(m-1))+lambda*gd(xi(i))+wij(i-1,m-1)+wij(i+1,m-1)+lambda*wij(i,m-2))/mu;
        if abs(wij(i,m-1)-z)>NORM
            NORM = abs(wij(i,m-1)-z);
        end
        wij(i,m-1) = z;
        %Step 9
        z=(-h^2*f(xi(n-1),yj(m-1))+gb(yj(m-1))+lambda*gd(xi(n-1))+wij(n-2,m-1)+lambda*wij(n-1,m-2))/mu;
        if abs (wij(n-1,m-1)-z) > NORM 
            NORM = abs (wij(n-1,m-1)-z);
        end
        wij(n-1,m-1)=z;

        %Step 10
        for j=m-2:-1:2
        
            %Step 11
            z=(-(h^2)*f(xi(1),yj(j))+ga(yj(j))+lambda*wij(1,j+1)+lambda*wij(1,j-1)+wij(2,j))/mu;
            if abs(wij(1,j)-z)>NORM
                NORM = abs(wij(1,j)-z);
            end
            wij(1,j)=z;

            %Step 12
            for i=2:n-2
                z=(-(h^2)*f(xi(i),yj(j))+wij(i-1,j)+lambda*wij(i,j+1)+wij(i+1,j)+lambda*wij(i,j-1))/mu;
                if abs (wij(i,j)-z)>NORM
                    NORM=abs(wij(i,j)-z);
                end
                wij(i,j)=z;
            end

            %STEP 13
            z=(-h^2*f(xi(n-1),yj(j))+gb(yj(j))+wij(n-2,j)+lambda*wij(n-1,j+1)+lambda*wij(n-1,j-1))/mu;
            if abs(wij(n-1,j)-z)>NORM
                NORM=abs(wij(n-1,j)-z);
            end
            wij(n-1,j)=z;
        end
    
        %Step 14
        z = (-h^2*f(xi(1),yj(1))+ga(yj(1))+lambda*gc(xi(1))+lambda*wij(1,2)+wij(2,1))/mu;
        if abs(wij(1,1)-z) > NORM
            NORM = abs(wij(1,1)-z);
        end
        wij(1,1) = z;

        %Step 15
        for i = 2:n-2
            z = (-h^2*f(xi(i),yj(1))+lambda*gc(xi(i))+wij(i-1,1)+lambda*wij(i,2)+wij(i+1,1))/mu;
            if abs(wij(i,1)-z)>NORM
                NORM = abs(wij(i,1)-z);
            end
            wij(i,1) = z;
        end

        %Step 16
        z=(-h^2*f(xi(n-1),yj(1))+gb(yj(1))+lambda*gc(xi(n-1))+wij(n-2,1)+lambda*wij(n-1,2))/mu;
            if abs (wij(n-1,1)-z) > NORM
                NORM= abs(wij(n-1,1)-z) ;
            end
            wij(n-1,1)=z;

        %Step 17
        if NORM<=TOL 
            %Step 18
            for i=1:n-1
                for j=1:m-1
                    disp("   "+double(xi(i))+"         "+double(yj(j))+"         "+double(wij(i,j)))
                end
            end
        %Step 19
        end
        %Step 20
        l=l+1;   
    end
end 

Answer 1

[Walter Roberson]

Yes, that is what you asked it to do. You display output whenever NORM <= TOL for any given iteration of l, and you iterate l up to N = 100, so you can have up to 100 displays of output.

Presumably something inside your loop is making changes so that the next iteration of l will be working with different inputs, so that the output would be different the next time... otherwise it would be pointless to continue iterating l once the tolerance had been met.