% FIGURE 4.7.1
%                              N = NUMBER OF POINTS IN X,Y DIRECTIONS
%                              A = PDE PARAMETER
      N = 40;
      A = 40.0;
      DX = 1.0/N;
      DY = DX;
%                              TRY VARIOUS VALUES OF THE PARAMETER OMEGA
      for W=1.0:0.05:1.95
%                              SET UNKNOWNS TO ZERO STARTING VALUES
         for I=2:N
         for J=2:N
            U(I,J) = 0.0;
         end
         end
%                              SET VALUES OF U ON BOUNDARY (THESE WILL
%                              NEVER CHANGE)
         for L=1:N+1
            X(L) = (L-1)*DX;
            Y(L) = (L-1)*DY;
            U(L,1) = 0.0;
            U(L,N+1) = X(L);
            U(1,L) = 0.0;
            U(N+1,L) = Y(L)^3;
         end
%                              BEGIN SOR ITERATION
         for ITER=1:10000
            ERMAX = 0.0;
%                              UPDATE UNKNOWNS USING SOR FORMULA
            for I=2:N
            for J=2:N
               U(I,J) = U(I,J) - 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  ...
               + A*U(I,J) - X(I)*Y(J)*(A*Y(J)^2-6) ) / (4/DX^2+A);
               ERR = abs(U(I,J) - X(I)*Y(J)^3);
               ERMAX = max(ERMAX,ERR);
            end
            end
            MAXITR = ITER;
%                              STOP IF ERROR < 10^(-6)
            if (ERMAX < 1.E-6)
               break
            end
         end
         W
         MAXITR
      end