4 nested for loops

I have a code with 4 nested for loops below. What i am trying to do is to calculate the error given by at 20th iteration for each and every omega value. As you can see from the code S depend on space(i and j), time(iteration) and omega(k) and thats why i used 4 nested for loops. When i try to output the error i just get 0 which is wrong since i created 4D matrix for S and i cannot display the summaries of variabes with more than 524288 elements i cannot figure out where it is going wrong.

What i want to ask is that if i am using for loops right for the scenario i explained above and if the loops are okay how can i display S to see where it goes wrong?

Any help you can give would be greatly appreciated.

clc
close all
clear all

N=21;
iteration=20;
deltax=1/(N-1);
deltay=deltax;
delomega=0.002;
omegai=1;
Nomega=1/(0.002);
x=linspace(0,1,21);
y=linspace(0,1,21);
omega=linspace(1,2,Nomega);

%analytical solution 
u=zeros(N,N,iteration);
for l=1:iteration
    for j=1:N
        for i=1:N
            u(i,j,l)=x(i)^2-y(j)^2;
        end
    end
end
s=zeros(N,N,iteration);
s(1,:,:)=u(1,:,:);
s(N,:,:)=u(N,:,:);
s(:,1,:)=u(:,1,:);
s(:,N,:)=u(:,N,:);

error=0;
for k=1:Nomega
    for l=1:iteration-1
        for j=2:N-1
            for i=2:N-1
                R(i,j,l+1)=s(i+1,j,l)+s(i-1,j,l+1)+s(i,j+1,l)+s(i,j-1,l+1)-4*(s(i,j,l))-(deltax^2)*s(i,j,l); 
                s(i,j,l+1,k)=s(i,j,l)+0.25*omega(k)*R(i,j,l+1);
            end  
        end
    end
     if l==iteration
         error=error+abs(u(i,j,l)-s(i,j,l));
     end
    error_(k)=error/((N-2)*(N-2));
end

plot(omega(:),error_(:));
xlabel('\omega');
ylabel('Error');
grid on
grid minor

[Stephen23]

I doubt that any human could make much sense of looking at an array with more than half a million elements... but it is trival to use indexing to look at any submatrix of your choice, and you could easily generate the indices from the data itself, e.g. find the first zero (or whatever the first "incorrect" data happens to be).