C***********************************************************************
C                                                                      *
C               POISSON FINITE-DIFFERENCE METHOD                       *
C                                                                      *
C***********************************************************************
C
C
C
C     TO APPROXIMATE THE SOLUTION TO THE POISSON EQUATION
C         DEL(U) = F(X,Y), A <= X <= B, C <= Y <= D,
C     SUBJECT TO BOUNDARY CONDITIONS:
C                     U(X,Y)=G(X,Y),
C     IF X = A OR X = B, FOR C <= Y <= D
C     IF Y = C OR Y = D, FOR A <= X <= B
C
C     INPUT:   ENDPOINTS A, B, C, D; INTEGERS M, N; TOLERENCE TOL;
C              MAXIMUM NUMBER OF ITERATIONS M.
C
C     OUTPUT:  APPROXIMATIONS W(I,J) TO U(X(I),Y(J)) FOR EACH
C              I=1,...,N-1 AND J=1,...,M-1 OR A MESSAGE THAT THE
C              MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.
C
C     INITIALIZATION
C     DIMENSION W(N-1,M-1),X(N-1),Y(M-1)
      DIMENSION W(5,4),X(5),Y(4)
      CHARACTER NAME*14,NAME1*14,AA*1
      INTEGER INP,OUP,FLAG,LBOUND
      LOGICAL OK
C     USE  -F(X,Y) INSTEAD OF F(X,Y) AND WE USED EXACT SOLUTION FOR
C     BOUNDARY VALUES
      F(XZ,YZ)=-XZ*EXP(YZ)
      G(XZ,YZ)= XZ*EXP(YZ)
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      WRITE(6,*) 'This is the Linear Finite Difference Method'
      WRITE(6,*) 'for Elliptic Equations.'
      WRITE(6,*) 'Have the functions F(x,y) and G(x,y) been created?'
      WRITE(6,*) 'Enter Y or N '
      WRITE(6,*) ' '
      READ(5,*)  AA
      IF(( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y' )) THEN
         OK = .FALSE.
19       IF (OK) GOTO 11
            WRITE(6,*) 'Input endpoints [A,B] on X-axis'
            WRITE(6,*) 'separated by blank.'
            WRITE(6,*) ' '
            READ(5,*) A, B
            WRITE(6,*) 'Input endpoints [C,D] on X-axis'
            WRITE(6,*) 'separated by blank.'
            WRITE(6,*) ' '
            READ(5,*) C, D
            IF ((A.GE.B).OR.(C.GE.D)) THEN
               WRITE(6,*) 'Left endpoint must be less'
               WRITE(6,*) 'than right endpoint'
            ELSE
               OK = .TRUE.
            ENDIF
         GOTO 19
11       OK = .FALSE.
12       IF (OK) GOTO 13
            WRITE(6,*) 'Input number of intervals n on the X-axis'
            WRITE(6,*) 'and m on the Y-axis separated by a blank.'
            WRITE(6,*) 'Note:  both n and m should be larger than 2.'
            WRITE(6,*) ' '
            READ(5,*) N, M
            IF ((N.LE.0).OR.(M.LE.0)) THEN
              WRITE(6,*) 'Must be positive integers '
            ELSE
              OK = .TRUE.
            ENDIF
         GOTO 12
13       OK = .FALSE.
14       IF (OK) GOTO 15
            WRITE(6,*) 'Input tolerance '
            WRITE(6,*) ' '
            READ(5,*) TOL
            IF (TOL.LE.0.0) THEN
               WRITE(6,*) 'Tolerance must be positive '
            ELSE
               OK = .TRUE.
            ENDIF
         GOTO 14
15       OK = .FALSE.
16       IF(OK) GOTO 17
            WRITE(6,*) 'Input the maximum number of iterations.'
            WRITE(6,*) ' '
            READ(5,*) LBOUND
            IF(LBOUND.LE.0) THEN
               WRITE(6,*) 'Must be a positive integer.'
            ELSE
               OK = .TRUE.
            ENDIF
            GO TO 16
17       CONTINUE
      ELSE
         WRITE(6,*) 'The program will end so that the functions'
         WRITE(6,*) 'F and G can be created.'
         OK = .FALSE.
      ENDIF
      IF(.NOT.OK) GOTO 400
      WRITE(6,*) 'Select output destinations: '
      WRITE(6,*) '1. Screen '
      WRITE(6,*) '2. Text file '
      WRITE(6,*) 'Enter 1 or 2 '
      WRITE(6,*) ' '
      READ(5,*) FLAG
      IF ( FLAG .EQ. 2 ) THEN
         WRITE(6,*) 'Input the file name in the form - '
         WRITE(6,*) 'drive:name.ext'
         WRITE(6,*) 'with the name contained within quotes'
         WRITE(6,*) 'as example:   ''A:OUTPUT.DTA'' '
         WRITE(6,*) ' '
         READ(5,*) NAME1
         OUP = 3
         OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW')
      ELSE
         OUP = 6
      ENDIF
      WRITE(OUP,*) 'POISSON EQUATION FINITE-DIFFERENCE METHOD'
      WRITE(OUP,*) A,B,C,D
      WRITE(OUP,*) TOL,M,N,LBOUND
      MM=M-1
      MMM=M-2
      NN=N-1
      NNN=N-2
C     XK WILL BE USED FOR K
C     STEP 1
      H=(B-A)/N
      XK=(D-C)/M
C     A,B,C,D WILL BE USED FOR X(0),X(N),Y(0),Y(M)
C     STEPS 2, 3 CONSTRUCT MESH POINTS
C     STEP 2
      DO 10 I=1,NN
10         X(I)=A+I*H
C     STEP 3
      DO 20 J=1,MM
20         Y(J)=C+J*XK
C     STEP 4
      DO 30 I=1,NN
           DO 30 J=1,MM
30              W(I,J)=0
C     STEP 5
C     USE V FOR LAMBDA, VV FOR MU
      V=H*H/(XK*XK)
      VV=2*(1+V)
      L=1
C     Z IS A NEW VALUE OF W(I,J) TO BE USED IN COMPUTING THE ERROR E
C     USE E INSTEAD OF NORM
C     STEP 6
100   IF (L.GT.LBOUND) GOTO 200
C     STEPS 7 THROUGH 20 PERFORM GAUSS-SEIDEL ITERATIONS
C          STEP 7
           Z=(H*H*F(X(1),Y(MM))+G(A,Y(MM))+V*G(X(1),D)+W(1,MM-1)*V+W(2,
     *     MM))/VV
           E=ABS(W(1,MM)-Z)
           W(1,MM)=Z
C          STEP 8
           DO 40 I=2,NNN
                Z=(H*H*F(X(I),Y(MM))+V*G(X(I),D)+W(I-1,MM)+W(I+1,MM)+V*
     *          W(I,MM-1))/VV
                IF(ABS(W(I,MM)-Z).GT.E) E=ABS(W(I,MM)-Z)
40         W(I,MM)=Z
C          STEP 9
           Z=(H*H*F(X(NN),Y(MM))+G(B,Y(MM))+V*G(X(NN),D)+W(NN-1,MM)+V*
     *     W(NN,MM-1))/VV
           IF(ABS(W(NN,MM)-Z).GT.E) E=ABS(W(NN,MM)-Z)
           W(NN,MM)=Z
C          STEP 10
           DO 50 LL=2,MMM
                J=MMM-LL+2
C               STEP 11
                Z=(H*H*F(X(1),Y(J))+G(A,Y(J))+V*W(1,J+1)+V*W(1,J-1)+W(2
     *          ,J))/VV
                IF(ABS(W(1,J)-Z).GT.E) E=ABS(W(1,J)-Z)
                W(1,J)=Z
C               STEP 12
                DO 60 I=2,NNN
                     Z=(H*H*F(X(I),Y(J))+W(I-1,J)+V*W(I,J+1)+V*W(I,J-1)
     *               +W(I+1,J))/VV
                     IF(ABS(W(I,J)-Z).GT.E) E=ABS(W(I,J)-Z)
60              W(I,J)=Z
C               STEP 13
                Z=(H*H*F(X(NN),Y(J))+G(B,Y(J))+W(NN-1,J)+V*W(NN,J+1)+V*
     *          W(NN,J-1))/VV
                IF(ABS(W(NN,J)-Z).GT.E) E=ABS(W(NN,J)-Z)
50         W(NN,J)=Z
C          STEP 14
           Z=(H*H*F(X(1),Y(1))+V*G(X(1),C)+G(A,Y(1))+V*W(1,2)+W(2,1))/VV
           IF(ABS(W(1,1)-Z).GT.E) E=ABS(W(1,1)-Z)
           W(1,1)=Z
C          STEP 15
           DO 70 I=2,NNN
                Z=(H*H*F(X(I),Y(1))+V*G(X(I),C)+W(I+1,1)+W(I-1,1)+V*W(I
     *          ,2))/VV
                IF(ABS(W(I,1)-Z).GT.E) E=ABS(W(I,1)-Z)
70         W(I,1)=Z
C          STEP 16
           Z=(H*H*F(X(NN),Y(1))+V*G(X(NN),C)+G(B,Y(1))+W(NN-1,1)+V*W(NN
     *     ,2))/VV
           IF(ABS(W(NN,1)-Z).GT.E) E=ABS(W(NN,1)-Z)
           W(NN,1)=Z
C          STEP 17
           IF(E.LE.TOL) THEN
                WRITE(OUP,2) L
C               STEP 18
                DO 80 I=1,NN
                DO 80 J=1,MM
80                   WRITE(OUP,3) I,J,X(I),Y(J),W(I,J)
C               STEP 19
                GOTO 400
           END IF
C          STEP 20
           L=L+1
      GOTO 100
C     STEP 21
200   WRITE(OUP,1) LBOUND
400   CLOSE(UNIT=5)
      CLOSE(UNIT=OUP)
      IF(OUP.NE.6) CLOSE(UNIT=6)
      STOP
1     FORMAT(1X,'METHOD FAILS AFTER ITERATION NO. ',I4)
2     FORMAT(1X,'NO. OF ITERATIONS IS',I4,3X,'ORDER OF OUTPUT IS I,J,X(I
     *),Y(J),W(I,J)')
3     FORMAT(1X,2(I3,2X),3(E15.8,3X))
      END
