% FIGURE 4.7.2
%                              N = NUMBER OF POINTS IN X,Y DIRECTIONS
      N = 40;
      DX = 2.0/N;
      DY = DX;
%                              TRY VARIOUS VALUES OF THE PARAMETER OMEGA
      for W=1.0:0.05:1.95
%                              SET UNKNOWNS AND BOUNDARY VALUES TO ZERO
         for I=1:N+1
         for J=1:N+1
            X(I) = -1.0 + (I-1)*DX;
            Y(J) = -1.0 + (J-1)*DY;
            U(I,J) = 0.0;
         end
         end
%                              BEGIN NLOR ITERATION
         for ITER=1:10000
            DELMAX = 0.0;
%                              UPDATE UNKNOWNS USING NLOR FORMULA
            for I=2:N
            LIM = 2;
            if (I <= N/2+1)
                LIM = N/2 + 2;
            end
            for J=LIM:N
               DELTA = W*((U(I+1,J)-2*U(I,J)+U(I-1,J))/DX^2 ...
                         +(U(I,J+1)-2*U(I,J)+U(I,J-1))/DY^2 ...
                                +(U(I+1,J)-U(I-1,J))/(2*DX) ...
                                +(U(I,J+1)-U(I,J-1))/(2*DY) ...
                        -X(I)^2*Y(J)^2*U(I,J)^3 - 1.0 ) ...
                  / (-2/DX^2 - 2/DY^2 - X(I)^2*Y(J)^2*3*U(I,J)^2);
               U(I,J) = U(I,J) - DELTA;
               DELMAX = max(DELMAX,abs(DELTA));
            end
            end
            MAXITR = ITER;
%                              STOP IF MAX CHANGE < 10^(-8)
            if (DELMAX < 1.E-8)
               break
            end
         end
         W
         MAXITR
      end