C
C PAGE 449-451: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
C
C FILE: SEIDEL.FOR
C
C ELLIPTIC PDE SOLVED BY DISCRETIZATION AND GAUSS-SEIDEL METHOD
C  (SEIDEL,F,G,BNDY,USTART,TRUE)
C
      PARAMETER  (IU = 9)   
      DIMENSION  U(IU,IU)   
      DATA  AX,BX/0.0,1.0/, AY,BY/0.0,1.0/  
      DATA  NX,NY/2*9/, ITMAX/20/     
      H = (BX - AX)/REAL(NX - 1)      
      DO 2 J = 1,NY 
        Y = AY + REAL(J-1)*H
        U(1,J)  = BNDY(AX,Y)
        U(NX,J) = BNDY(BX,Y)
   2  CONTINUE
      DO 3 I = 1,NX 
        X = AX + REAL(I-1)*H
        U(I,1)  = BNDY(X,AY)
        U(I,NY) = BNDY(X,BY)
   3  CONTINUE
      DO 5 J = 2,NY-1       
        Y = AY + REAL(J-1)*H
        DO 4 I = 2,NX-1     
          X = AX + REAL(I-1)*H
          U(I,J) = USTART(X,Y)
   4    CONTINUE
   5  CONTINUE    
      PRINT 8,0,((U(I,J),I = 1,NX),J = 1,NY)    
      CALL SEIDEL(AX,AY,NX,NY,H,ITMAX,U,IU)     
      PRINT 8,ITMAX,((U(I,J),I = 1,NX),J = 1,NY)
      DO 7 J = 1,NY 
        Y = AY + REAL(J-1)*H
        DO 6 I = 1,NX       
          X = AX + REAL(I-1)*H
          U(I,J) = ABS( TRUE(X,Y) - U(I,J) )  
   6    CONTINUE
   7  CONTINUE    
      PRINT 8,ITMAX,((U(I,J),I = 1,NX),J = 1,NY)
   8  FORMAT(//4X,I5,//(9(2X,E12.5))) 
      STOP
      END 
        
      FUNCTION F(X,Y)       
      F = -25.0 
      RETURN      
      END 
        
      FUNCTION G(X,Y)       
      G = 0.0 
      RETURN      
      END 
        
      FUNCTION BNDY(X,Y)    
      BNDY = TRUE(X,Y)    
      RETURN      
      END 
        
      FUNCTION USTART(X,Y)  
      USTART = 1.0
      RETURN      
      END 
        
      FUNCTION TRUE(X,Y)    
      TRUE = 0.5*(COSH(5.0*X) + COSH(5.0*Y))/COSH(5.0)
      RETURN      
      END 
  
      SUBROUTINE SEIDEL(AX,AY,NX,NY,H,ITMAX,U,IU) 
      DIMENSION  U(IU,NY)   
      HSQ = H*H   
      DO 4 K = 1,ITMAX      
        DO 3 J = 2,NY-1     
          Y = AY + REAL(J-1)*H
          DO 2 I = 2,NX-1   
            X = AX + REAL(I-1)*H      
            V = U(I+1,J) + U(I-1,J) + U(I,J+1) + U(I,J-1) 
            U(I,J) = (V - HSQ*G(X,Y))/(4.0 - HSQ*F(X,Y))  
   2      CONTINUE
   3    CONTINUE  
   4  CONTINUE    
      RETURN      
      END