*DECK PRO4P9
      SUBROUTINE INITAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
     5     /FE/  AL(N),  BM(N),  CU(N), BRHS(N)
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  LENGTH
      L=1.0D0
C...
C...  GRID SPACING
      DX=L/DFLOAT(N-1)
C...
C...  SET UP THE COEFFICIENT MATRIX IN BAND STORAGE MODE. THIS IS
C...  DONE ONLY ONCE SINCE THE COEFFICIENT MATRIX IS CONSTANT
C...
C...  LOWER DIAGONAL
      DO 2 I=1,N
         IF(I.EQ.1)THEN
            AL(1)=0.0D0
         ELSE
     +   IF(I.EQ.N)THEN
            AL(N)=1.0D0/6.0D0
         ELSE
            AL(I)=1.0D0/6.0D0
         END IF
2     CONTINUE
C...
C...  MAIN DIAGONAL
      DO 3 I=1,N
         IF(I.EQ.1)THEN
            BM(1)=2.0D0/6.0D0
         ELSE
     +   IF(I.EQ.N)THEN
            BM(N)=2.0D0/6.0D0
         ELSE
            BM(I)=4.0D0/6.0D0
         END IF
3     CONTINUE
C...
C...  UPPER DIAGONAL
      DO 4 I=1,N
         IF(I.EQ.1)THEN
            CU(1)=1.0D0/6.0D0
         ELSE
     +   IF(I.EQ.N)THEN
            CU(N)=0.0D0
         ELSE
            CU(I)=1.0D0/6.0D0
         END IF
4     CONTINUE
C...
C...  INITIAL CONDITION
      DO 1 I=1,N
         X(I)=DFLOAT(I-1)*DX
         IF(I.LT.(N+1)/2)U(I)=0.0D0
         IF(I.GT.(N+1)/2)U(I)=1.0D0
         IF(I.EQ.(N+1)/2)U(I)=0.5D0
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
     5     /FE/  AL(N),  BM(N),  CU(N), BRHS(N)
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  LEFT BOUNDARY CONDITION
      H1=0.0D0
C...
C...  RIGHT BOUNDARY CONDITION
      HN=0.0D0
C...
C...  RIGHT HAND SIDE VECTOR (WITH D = 1)
      DO 1 I=1,N
         IF(I.EQ.1)THEN
            BRHS(I)=-H1/DX+(U(I+1)-U(I))/DX**2
         ELSE
     +   IF(I.EQ.N)THEN
            BRHS(I)= HN/DX+(U(I-1)-U(I))/DX**2
         ELSE
            BRHS(I)=(U(I+1)-2.0D0*U(I)+U(I-1))/DX**2
         END IF
1     CONTINUE
C...
C...  SOLVE THE LINEAR ALGEBRAIC EQUATIONS BY SUBROUTINE TRIDAG, WHICH
C...  RETURNS THE DERIVATIVE VECTOR UT IN COMMON/F/
      CALL TRIDAG(AL,BM,CU,BRHS,UT,N)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
     5     /FE/  AL(N),  BM(N),  CU(N), BRHS(N)
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DIMENSION ARRAYS FOR THE ANALYTICAL SOLUTION
      DIMENSION UE(N), DIFF(N)
C...
C...  PRINT THE COEFFICIENT MATRIX FOR VERIFICATION
      IP=IP+1
      IF(IP.EQ.1)THEN
         WRITE(NO,13)
13       FORMAT(/,' COEFFICIENT MATRIX',/)
         DO 14 I=1,N
            WRITE(NO,15)I,AL(I),BM(I),CU(I)
15          FORMAT(I5,3F12.4)
14       CONTINUE
         WRITE(NO,16)
16       FORMAT(//)
      END IF
C...
C...  COMPUTE THE TOTAL ENERGY (RHO = CP = 1)
      CALL DERV
      ENERGY=0.0D0
      DO 1 I=2,N
         ENERGY=ENERGY+0.5D0*(U(I)+U(I-1))*DX
1     CONTINUE
C...
C...  CALCULATE THE EXACT SOLUTION, AND THE DIFFERENCE BETWEEN THE
C...  NUMERICAL AND EXACT SOLUTIONS
      DO 3 I=1,N,10
         UE(I)=SERIES(I)
         DIFF(I)=U(I)-UE(I)
3     CONTINUE
C...
C...  PRINT THE NUMERICAL AND EXACT SOLUTIONS TO EQUATIONS (1) TO (4)
      WRITE(NO,2)T,ENERGY,(U(I),I=1,N,10),
     1                   (UE(I),I=1,N,10),
     2                 (DIFF(I),I=1,N,10)
2     FORMAT('    T = ',F6.2,'  ENERGY = ',F9.5,
     1           /,14X,'   X=0',5X,'X=0.2',5X,'X=0.4',
     2               5X,'X=0.6',5X,'X=0.8',5X,'X=1.0'/,
     3                         '    U(X,T)',6F10.6,/,
     4                         '   UE(X,T)',6F10.6,/,
     5                         ' DIFF(X,T)',6F10.6,/)
      RETURN
      END
      DOUBLE PRECISION FUNCTION SERIES(IG)
C...
C...  FUNCTION SERIES EVALUATES THE ANALYTICAL SOLUTION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=51)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/   U(N)
     2      /F/  UT(N)
     3     /XG/      L,     DX,   X(N)
     4      /I/     IP
C...
C...  TYPE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION L
C...
C...  DEFINE MAXIMUM NUMBER OF TERMS IN THE SERIES (THIS IS SET TO A
C...  HIGH VALUE TO ACHIEVE CONVERGENCE TO THREE FIGURES FOR T = 0;
C...  BEYOND T = 0, THE SERIES CONVERGES RAPIDLY)
      PARAMETER (IS=1000)
C...
C...  PI
      PI=4.0D0*DATAN(1.0D0)
C...
C...  EVALUATE SERIES SOLUTION
      SUM=0.25D0
      SIGN=1.0D0
      DO 1 I=1,IS,2
         EI=DFLOAT(I)*PI
         EXI=-(EI**2)*T
         SIGN=-1.0D0*SIGN
         CI=SIGN/EI
         SUM=SUM+CI*DEXP(EXI)*DCOS(EI*X(IG))
C...
C...     AVOID UNDERFLOW OF THE EXP FUNCTION
         IF(EXI.LT.-100.0D0)GO TO 2
1     CONTINUE
C...
C...  SUMMATION OF SERIES IS COMPLETE
2     SERIES=2.0*SUM
      END
