C FIGURE 4.8.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),R(0:N,0:N),P(0:N,0:N),                       
     &         AP(0:N,0:N),X(0:N),Y(0:N)
      DOUBLE PRECISION LAMBDA                                           
      DX = 1.D0/N                                                       
      DY = DX                                                           
C                              U = 0                                    
      DO 5 I=1,N-1                                                      
      DO 5 J=1,N-1                                                      
         U(I,J) = 0.0                                                      
    5 CONTINUE
C                              SET VALUES OF U,P 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                                                   
         P(L,0) = 0.0                                                   
         U(L,N) = X(L)                                                  
         P(L,N) = 0.0                                                  
         U(0,L) = 0.0                                                   
         P(0,L) = 0.0                                                   
         U(N,L) = Y(L)**3                                               
         P(N,L) = 0.0                                               
   10 CONTINUE                                                          
C                              R = B - A*U
C                              P = R                                    
      DO 15 I=1,N-1                                                     
      DO 15 J=1,N-1                                                     
         R(I,J) = (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)
         P(I,J) = R(I,J)                                                
   15 CONTINUE                                                          
C                              BEGIN CONJUGATE GRADIENT ITERATION       
      DO 30 ITER=1,1000                                                 
C                              AP = A*P                         
         DO 20 I=1,N-1                                                  
         DO 20 J=1,N-1                                                  
            AP(I,J) = (-P(I+1,J)+2*P(I,J)-P(I-1,J))/DX**2    
     &               +(-P(I,J+1)+2*P(I,J)-P(I,J-1))/DY**2    
     &               + A*P(I,J)   
   20    CONTINUE                                                       
C                              LAMBDA = R*P / P*AP                      
         LAMBDA = DOT(R,P,N)/DOT(P,AP,N)                                
C                              U = U + LAMBDA*P                         
         CALL UPDATE(U,U,LAMBDA,P,N)                                    
C                              R = R - LAMBDA*AP                        
         CALL UPDATE(R,R,-LAMBDA,AP,N)                                  
C                              ALPHA = -R*AP / P*AP                     
         ALPHA = -DOT(R,AP,N)/DOT(P,AP,N)                               
C                              P = R + ALPHA*P                          
         CALL UPDATE(P,R,ALPHA,P,N)                                     
C                              CALCULATE MAXIMUM ERROR                  
         ERMAX = 0.0                                                    
         DO 25 I=1,N-1                                                  
         DO 25 J=1,N-1                                                  
            ERR = ABS(U(I,J) - X(I)*Y(J)**3)                            
            ERMAX = MAX(ERMAX,ERR)                                      
   25    CONTINUE                                                       
         MAXITR = ITER                                                  
C                              STOP IF ERROR < 10**(-6)              
         IF (ERMAX.LT.1.D-6) GO TO 35                                   
   30 CONTINUE                                                          
   35 PRINT 40, MAXITR                                                  
   40 FORMAT (I10)                                                      
      STOP                                                              
      END                                                               
      FUNCTION DOT(X,Y,N)                                               
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION X(0:N,0:N),Y(0:N,0:N)                                   
C                              CALCULATE DOT PRODUCT OF X AND Y         
      DOT = 0.0                                                         
      DO 5 I=1,N-1                                                      
      DO 5 J=1,N-1                                                      
         DOT = DOT + X(I,J)*Y(I,J)                                         
    5 CONTINUE
      RETURN                                                            
      END                                                               
      SUBROUTINE UPDATE(X,Y,S,Z,N)                                      
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      DIMENSION X(0:N,0:N),Y(0:N,0:N),Z(0:N,0:N)                        
C                              CALCULATE X = Y + S*Z                    
      DO 5 I=1,N-1                                                      
      DO 5 J=1,N-1                                                      
         X(I,J) = Y(I,J) + S*Z(I,J)                                        
    5 CONTINUE
      RETURN                                                            
      END