*DECK PRO2P7
      SUBROUTINE INITAL
C...
C...  SUBROUTINE INITAL IS CALLED BY MAIN PROGRAM PRO1P3 TO DEFINE THE
C...  INITIAL CONDITION IN THE NUMOL SOLUTION OF EQUATIONS (2.37) TO
C...  (2.43) (VIA EQUATIONS (2.45) TO (2.47))
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(21,21)
     2      /F/ UT(21,21)
     3      /S/UXX(21,21),UYY(21,21),     X(21),     Y(21),
     4             US(21),   USX(21),   USY(21),  USXX(21),  USYY(21),
     5                 PI,      C(3),        XL,        YL
     6      /I/        NX,        NY
C...
C...  CONSTANTS USED IN THE THREE PDE PROBLEMS
      PI=ACOS(-1.)
      C(1)=PI
      C(2)=SQRT(PI**2+1.0)
      C(3)=SQRT(PI**2-1.0)
C...
C...  SPATIAL PARAMETERS
      NX=21
      NY=21
      XL=1.0
      YL=1.0
      DX=XL/FLOAT(NX-1)
      DY=YL/FLOAT(NY-1)
C...
C...  INITIALIZATION OVER THE SPATIAL GRID
      DO 1 I=1,NX
      DO 1 J=1,NY
C...
C...     SPATIAL GRID
         X(I)=DX*FLOAT(I-1)
         Y(J)=DY*FLOAT(J-1)
C...
C...     INITIAL CONDITION (2.48)
         U(I,J)=1.0
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV
C...
C...  SUBROUTINE DERV IS CALLED BY INTEGRATOR RKF45 TO DEFINE THE
C...  ODES IN THE NUMOL SOLUTION OF EQUATIONS (2.45) TO (2.47)
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(21,21)
     2      /F/ UT(21,21)
     3      /S/UXX(21,21),UYY(21,21),     X(21),     Y(21),
     4             US(21),   USX(21),   USY(21),  USXX(21),  USYY(21),
     5                 PI,      C(3),        XL,        YL
     6      /I/        NX,        NY
C...
C...  BOUNDARY CONDITION AT X = 0, EQUATION (2.40)
      DO 1 J=1,NY
1     U(1,J)=0.
C...
C...  BOUNDARY CONDITION AT X = 1, EQUATION (2.41)
      DO 2 J=1,NY
2     U(NX,J)=SINH(C(NORUN)*X(NX))*SIN(PI*Y(J))
C...
C...  BOUNDARY CONDITION AT Y = 0, EQUATION (2.42)
      DO 3 I=1,NX
3     U(I,1)=0.
C...
C...  BOUNDARY CONDITION AT Y = 1, EQUATION (2.43)
      DO 4 I=1,NX
4     U(I,NY)=0.
C...
C...  UXX
C...
C...     0 LE Y LE 1
         DO 5 J=1,NY
C...
C...     0 LE X LE 1
         DO 6 I=1,NX
C...
C...     TRANSFER U(I,J) TO US(I)
6        US(I)=U(I,J)
C...
C...        USX
            CALL DSS002(0.,XL,NX,US,USX)
C...
C...        USXX
            CALL DSS002(0.,XL,NX,USX,USXX)
C...
C...     TRANSFER USXX TO UXX
         DO 7 I=1,NX
7        UXX(I,J)=USXX(I)
5     CONTINUE
C...
C...  UYY
C...
C...     0 LE X LE 1
         DO 8 I=1,NX
C...
C...     0 LE Y LE 1
         DO 9 J=1,NY
C...
C...     TRANSFER U(I,J) TO US(I)
9        US(J)=U(I,J)
C...
C...        USY
            CALL DSS002(0.,YL,NY,US,USY)
C...
C...        USYY
            CALL DSS002(0.,YL,NY,USY,USYY)
C...
C...     TRANSFER USYY TO UYY
         DO 10 J=1,NY
10       UYY(I,J)=USYY(J)
8        CONTINUE
C...
C...  LAPLACE'S EQUATION VIA EQUATION (2.45)
      IF(NORUN.EQ.1)THEN
         DO 11 I=1,NX
         DO 11 J=1,NY
            UT(I,J)=UXX(I,J)+UYY(I,J)
11       CONTINUE
C...
C...  POISSON'S EQUATION VIA EQUATION (2.46)
      ELSE IF(NORUN.EQ.2)THEN
         DO 12 I=1,NX
         DO 12 J=1,NY
            UT(I,J)=UXX(I,J)+UYY(I,J)-SINH(C(NORUN)*X(I))*SIN(PI*Y(J))
 12      CONTINUE
C...
C...  HELMHOLTZ'S EQUATION VIA EQUATION (2.47)
      ELSE IF(NORUN.EQ.3)THEN
         DO 13 I=1,NX
         DO 13 J=1,NY
            UT(I,J)=UXX(I,J)+UYY(I,J)+U(I,J)
13       CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  SUBROUTINE PRINT IS CALLED BY MAIN PROGRAM PRO1P3 TO PRINT THE
C...  NUMOL SOLUTION OF EQUATIONS (2.37) TO (2.43)
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(21,21)
     2      /F/ UT(21,21)
     3      /S/UXX(21,21),UYY(21,21),     X(21),     Y(21),
     4             US(21),   USX(21),   USY(21),  USXX(21),  USYY(21),
     5                 PI,      C(3),        XL,        YL
     6      /I/        NX,        NY
C...
C...  DIMENSION THE ARRAYS FOR THE EXACT SOLUTION
      REAL      UE(21,21),DIFF(21,21)
C...
C...  CALCULATE THE EXACT SOLUTION, AND THE DIFFERENCE BETWEEN THE
C...  NUMOL AND EXACT SOLUTIONS
      DO 4 J=NY,1,-4
      DO 4 I=1,NX,4
      UE(I,J)=EXACT(I,J)
      DIFF(I,J)=U(I,J)-UE(I,J)
4     CONTINUE
C...
C...  PRINT THE NUMOL SOLUTION TO EQUATIONS (2.45) TO (2.47)
      WRITE(NO,1)T
1     FORMAT(' "TIME" = ',F6.2,/,20X,'  X=0',4X,'X=0.2',4X,'X=0.4',4X,
     1                             'X=0.6',4X,'X=0.8',4X,'  X=1',/)
      DO 3 J=NY,1,-4
      WRITE(NO,2)Y(J),( U(I,J),I=1,NX,4),
     1                (UE(I,J),I=1,NX,4),
     2              (DIFF(I,J),I=1,NX,4),
     3                (UT(I,J),I=1,NX,4)
2     FORMAT(' Y = ',F3.1,'  U(X,Y)',6F9.5,/,
     1  9X,'UE(X,Y)',6F9.5,/,12X,'DIFF',6F9.5,/,
     2  9X,'UT(X,Y)',6F9.2,/)
3     CONTINUE
      RETURN
      END
      REAL FUNCTION EXACT(I,J)
C...
C...  FUNCTION EXACT COMPUTES THE EXACT SOLUTION TO EQUATIONS (2.37)
C...  TO (2.43), I.E., EQUATION (2.44), AT GRID INDICES I AND J
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(21,21)
     2      /F/ UT(21,21)
     3      /S/UXX(21,21),UYY(21,21),     X(21),     Y(21),
     4             US(21),   USX(21),   USY(21),  USXX(21),  USYY(21),
     5                 PI,      C(3),        XL,        YL
     6      /I/        NX,        NY
C...
C...  EXACT SOLUTION AT X AND Y
      EXACT=SINH(C(NORUN)*X(I))*SIN(PI*Y(J))
      RETURN
      END
