C FIGURE 4.7.1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              N = NUMBER OF POINTS IN X,Y DIRECTIONS   
C                              A = PDE PARAMETER                        
      PARAMETER (N=40, A=40.0)                                          
      DIMENSION U(0:N,0:N),X(0:N),Y(0:N)                                
      DX = 1.D0/N                                                       
      DY = DX                                                           
C                              TRY VARIOUS VALUES OF THE PARAMETER OMEGA
      DO 35 IW=0,19                                                     
         W = 1.0 + IW*0.05                                              
C                              SET UNKNOWNS TO ZERO STARTING VALUES     
         DO 5 I=1,N-1                                                   
         DO 5 J=1,N-1                                                   
            U(I,J) = 0.0                                                   
    5    CONTINUE
C                              SET VALUES OF U ON BOUNDARY (THESE WILL  
C                              NEVER CHANGE)                            
         DO 10 L=0,N                                                    
            X(L) = L*DX                                                 
            Y(L) = L*DY                                                 
            U(L,0) = 0.0                                                
            U(L,N) = X(L)                                               
            U(0,L) = 0.0                                                
            U(N,L) = Y(L)**3                                            
   10    CONTINUE                                                       
C                              BEGIN SOR ITERATION                      
         DO 20 ITER=1,10000                                             
            ERMAX = 0.0                                                 
C                              UPDATE UNKNOWNS USING SOR FORMULA        
            DO 15 I=1,N-1                                               
            DO 15 J=1,N-1                                               
               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) ) / (2/DX**2+2/DY**2+A)     
               ERR = ABS(U(I,J) - X(I)*Y(J)**3)                         
               ERMAX = MAX(ERMAX,ERR)                                   
   15       CONTINUE                                                    
            MAXITR = ITER                                               
C                              STOP IF ERROR < 10**(-6)              
            IF (ERMAX.LT.1.D-6) GO TO 25                                
   20    CONTINUE                                                       
   25    PRINT 30, W,MAXITR                                             
   30    FORMAT (5X,F8.2,I10)                                           
   35 CONTINUE                                                          
      STOP                                                              
      END