C FIGURE 4.7.2
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              N = NUMBER OF POINTS IN X,Y DIRECTIONS   
      PARAMETER (N=40)                                                  
      DIMENSION U(0:N,0:N),X(0:N),Y(0:N)                                
      DX = 2.D0/N                                                       
      DY = DX                                                           
C                              TRY VARIOUS VALUES OF THE PARAMETER OMEGA
      DO 30 IW=0,19                                                     
         W = 1.0 + IW*0.05                                              
C                              SET UNKNOWNS AND BOUNDARY VALUES TO ZERO 
         DO 5 I=0,N                                                     
         DO 5 J=0,N                                                     
            X(I) = -1.0 + I*DX                                          
            Y(J) = -1.0 + J*DY                                          
            U(I,J) = 0.0                                                   
    5    CONTINUE
C                              BEGIN NLOR ITERATION                     
         DO 15 ITER=1,10000                                             
            DELMAX = 0.0                                                
C                              UPDATE UNKNOWNS USING NLOR FORMULA       
            DO 10 I=1,N-1                                               
            LIM = 1                                                     
            IF (I.LE.N/2) LIM = N/2 + 1                                 
            DO 10 J=LIM,N-1                                             
               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))                          
   10       CONTINUE                                                    
            MAXITR = ITER                                               
C                              STOP IF MAX CHANGE < 10**(-8)         
            IF (DELMAX.LT.1.D-8) GO TO 20                               
   15    CONTINUE                                                       
   20    PRINT 25, W,MAXITR                                             
   25    FORMAT (5X,F8.2,I10)                                           
   30 CONTINUE                                                          
      STOP                                                              
      END