      SUBROUTINE INITAL
C...
C...  ONE-DIMENSIONAL CONVECTION AND DIFFUSION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  SPATIAL DIFFERENTIATOR
C...
C...     NCASE = 1 - DSS020
C...
C...     NCASE = 2 - DSS024
C...
      NCASE=1
C...
C...  SYSTEM LENGTH
      ZL=1.0D0
C...
C...  SYSTEM PARAMETERS
      PE=10.0D0
C...
C...  INITIAL CONDITIONS
      DO 1 I=1,NZ
         C(I)=0.0D0
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
C...
C...  INITIALIZE COUNTERS FOR PRINTING AND PLOTTING
      IP=0
      NP=81
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  CZ
      IF(NCASE.EQ.1)CALL DSS020(0.0D0,ZL,NZ,C,CZ,1.0D0)
      IF(NCASE.EQ.2)CALL DSS024(0.0D0,ZL,NZ,C,CZ,1.0D0)
C...
C...  BOUNDARY CONDITION, Z = 0
         NL=2
         CZ(1)=PE*(C(1)-FB(T))
C...
C...  BOUNDARY CONDITION, Z = 1
         NU=2
         CZ(NZ)=0.0D0
C...
C...  CZZ
      CALL DSS044(0.0D0,ZL,NZ,C,CZ,CZZ,NL,NU)
C...
C...  PDE
      DO 1 I=1,NZ
         CT(I)=CZZ(I)-PE*CZ(I)
 1    CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION FB(T)
C...
C...  UNIT STEP FUNCTION
      IF(T.LT.0.0D0)THEN
         FB=0.0D0
      ELSE
     +IF(T.GE.0.0D0)THEN
         FB=1.0D0
      END IF
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  MONITOR THE OUTPUT
      WRITE(*,*)IP
C...
C...  WRITE NUMERICAL AND SERIES SOLUTIONS EVERY 40TH CALL TO PRINT
      IP=IP+1
      IF (((IP-1)/40*40).EQ.(IP-1))THEN
         WRITE(NO,1)T,(C(I),I=1,NZ)
1        FORMAT(/,' t = ',F5.2,/,F10.5,/,(5F10.5))
      END IF
C...
C...  MATLAB PLOTTING OF THE SOLUTION
      CALL PLOTM
C...
C...  SERIES SOLUTION AT END OF NUMERICAL SOLUTION
      IF(IP.EQ.NP)CALL SERIES(NI,NO)
      RETURN
      END
      SUBROUTINE PLOTM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  DEFINE A FILE FOR MATLAB PLOTTING OF THE NUMERICAL SOLUTION
      IF(IP.EQ.1)OPEN(1,FILE='ma2b.mat')
C...
C...  WRITE SOLUTION AT Z = 1 FOR MATLAB PLOTTING (EVERY FIFTH POINT IS
C...  USED SO THAT THE INDIVIDUAL PLOTTED POINTS CAN BE DISTINGUISHED)
      IF(((IP-1)/5*5).EQ.(IP-1))WRITE(1,1)T,C(NZ)
1     FORMAT(2F10.4)
      RETURN
      END
      SUBROUTINE SERIES(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  DEFINE A FILE FOR MATLAB PLOTTING OF THE ANALYTICAL SOLUTION
      OPEN(2,FILE='ma2a.out')
C...
C...  NUMBER OF EIGENVALUES
      NS=10
C...
C...  NUMBER OF PANELS IN THE NUMERICAL QUADRATURE (INTEGRATION)
      NQ=1000
C...
C...  LIMITS OF INTEGRATION
      ZL=0.0D0
      ZU=1.0D0
C...
C...  SYSTEM LENGTH
      ZLU=1.0D0
C...
C...  SYSTEM PARAMETERS
      PE=10.0D0
C...
C...  COMPUTE AND PRINT NS EIGENVALUES
      CALL EIGVAL
      WRITE(NO,1)
      WRITE( *,1)
1     FORMAT(//,' Eigenvalues',/)
      WRITE(NO,2)(I,ZMU(I),I=1,NS)
2     FORMAT(I5,F10.5)
C...
C...  TEST ORTHOGONALITY OF EIGENFUNCTIONS BY COMPUTING THE INTEGRAL
C...  OF THE PRODUCT OF TWO EIGENFUNCTIONS FOR I AND J = 1, 2,..., NS
      WRITE(NO,13)
      WRITE( *,13)
13    FORMAT(//,' Product integrals',/)
      DO 10 I=1,NS
      DO 11 J=1,NS
         CALL  ORTHO(I,J,PIJ)
         WRITE(NO,12)I,J,PIJ
12       FORMAT(' i = ',I3,'   j = ',I3,'   Pij = ',F10.5)
         IF(J.EQ.NS)WRITE(NO,14)
14       FORMAT(/)
11    CONTINUE
10    CONTINUE
C...
C...  COMPUTE AND PRINT NS FOURIER COEFFICIENTS
      NS=200
      CALL EIGVAL
      WRITE(NO,22)
      WRITE( *,22)
22    FORMAT(//,' Fourier coefficients',/)
      DO 20 I=1,NS
         CALL COEFF(I)
         WRITE(NO,21)I,CF(I)
21       FORMAT(I5,F10.7)
20       CONTINUE
C...
C...  FOURIER SERIES FOR UNITY WITH NS TERMS AT A SERIES OF VALUES OF Z
      NS=200
      WRITE(NO,30)
      WRITE( *,30)
30    FORMAT(//,' Fourier series for unity',/)
C...
C...  STEP THROUGH NZV VALUES OF Z
      NZV=11
      DO 32 J=1,NZV
C...
C...     SET Z
         Z=DFLOAT(J-1)/DFLOAT(NZV-1)
C...
C...     COMPUTE AND PRINT SERIES FOR UNITY
         ONE=UNITY(Z)
         WRITE(NO,31)Z,ONE
31       FORMAT(F10.2,F10.5)
32    CONTINUE
C...
C...  EVALUATE AND PRINT THE SERIES SOLUTION AT Z = 1 FOR A SERIES
C...  OF VALUES OF T USING NS TERMS IN THE SERIES
      WRITE(NO,42)
      WRITE( *,42)
42    FORMAT(//,' Series solution at z = 1, phi1(1,t)',/)
      Z=1.0D0
      DT=0.0025D0
      NTM=81
      DO 40 I=1,NTM
C...
C...     T AND SERIES SOLUTION
         T=DT*DFLOAT(I-1)
         UZT=PHI1(Z,T)
C...
C...     PRINT SOLUTION
         WRITE(NO,41)Z,T,UZT
41       FORMAT(F10.3,F10.4,F10.5)
C...
C...     WRITE SOLUTION FOR MATLAB PLOTTING
         WRITE(2,43)T,UZT
43       FORMAT(2F10.4)
40    CONTINUE
      RETURN
      END
      SUBROUTINE EIGVAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  FUNCTION FOR THE TRANSCENDENTAL EQUATION THAT DEFINES THE
C...  EIGENVALUES
      EXTERNAL F1
C...
C...  PI
      PI=4.0*DATAN(1.0D0)
C...
C...  ERROR TOLERANCE
      EPS=1.0D-06
C...
C...  COMPUTE NS EIGENVALUES
      DO 1 I=1,NS
C...
C...     BOUND ROOT (EIGENVALUE)
         ZP=DFLOAT(I-1)*PI
         ZN=DFLOAT(I  )*PI
C...
C...     COMPUTE ROOT AND STORE AS AN EIGENVALUE
         ZMU(I)=ROOT(F1,ZN,ZP,EPS)
C...
C...  NEXT EIGENVALUE
1     CONTINUE
      RETURN
      END
      SUBROUTINE ORTHO(I,J,PIJ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  EXTERNAL THE FUNCTION FOR THE INTEGRAND OF THE PRODUCT INTEGRAL
      EXTERNAL F2
C...
C...  SET THE INDICES OF THE EIGENFUNCTIONS
      IN=I
      JN=J
C...
C...  PERFORM THE INTEGRATION OF THE PRODUCT OF TWO EIGENFUNCTIONS
C...  FROM ZL = 0 TO ZU = 1
      PIJ=SIMP(F2,ZL,ZU,NQ)
      RETURN
      END
      SUBROUTINE COEFF(I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  EXTERNAL THE FUNCTION FOR THE INTEGRAND OF THE NUMERATOR AND
C...  DENOMINATOR INTEGRALS OF THE FOURIER COEFFICIENT
      EXTERNAL F3, F4
C...
C...  SET THE INDEX OF THE FOURIER COEFFICIENT
      IN=I
C...
C...  PERFORM THE INTEGRATION FOR THE NUMERATOR AND DENOMINATOR
C...  INTEGRALS FROM ZL = 0 TO ZU = 1
      ZNUM=SIMP(F3,ZL,ZU,NQ)
      ZDEN=SIMP(F4,ZL,ZU,NQ)
C...
C...  FOURIER COEFFICIENT
      CF(I)=ZNUM/ZDEN
      RETURN
      END
      DOUBLE PRECISION FUNCTION UNITY(Z)
C...
C...  FUNCTION UNITY COMPUTES THE SERIES EXPANSION OF UNITY (THE
C...  NUMBER ONE) AT Z
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  INITIALIZE SERIES
      SUM=0.0D0
C...
C...  PRECOMPUTE THE COMMON EXPONENTIAL FACTOR
      EXPZ=DEXP((PE/2.0D0)*Z)
C...
C...  COMPUTE SERIES WITH NS TERMS
      DO 1 I=1,NS
C...
C...     EIGENFUNCTION PSI(I)
         PSII=PE/(2.0D0*ZMU(I))*DSIN(ZMU(I)*Z)+DCOS(ZMU(I)*Z)
C...
C...     RUNNING SUM
         SUM=SUM+CF(I)*PSII
C...
C...  CONTINUE SUMMATION
1     CONTINUE
C...
C...  INCLUDE EXPONENTIAL FACTOR
      UNITY=EXPZ*SUM
C...
C...  CALCULATION IS COMPLETE
      RETURN
      END
      DOUBLE PRECISION FUNCTION PHI1(Z,T)
C...
C...  FUNCTION PHI1 COMPUTES THE SERIES SOLUTION AT Z AND T
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/      TM,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  INITIALIZE SERIES
      SUM=0.0D0
C...
C...  PRECOMPUTE THE COMMON EXPONENTIAL FACTOR
      EXPZ=DEXP((PE/2.0D0)*Z)
C...
C...  COMPUTE SERIES WITH NS TERMS
      DO 1 I=1,NS
C...
C...     EXPONENTIAL IN T
         TLAM=ZMU(I)**2+PE**2/4.0D0
C...
C...     TEST MAGNITUDE OF TLAM*T TO AVOID UNDERFLOW OF EXP FUNCTION
         TLAM=-TLAM*T
         IF(TLAM.LT.-1.0D+20)THEN
C...
C...     TERMS ARE NEGLIGIBLY SMALL SO TERMINATE SUMMING SERIES
            GO TO 2
         ELSE
C...
C...     CONTINUE TO SUM SERIES
            EXPT=DEXP(TLAM)
         END IF
C...
C...     EIGENFUNCTION PSI(I)
         PSII=PE/(2.0D0*ZMU(I))*DSIN(ZMU(I)*Z)+DCOS(ZMU(I)*Z)
C...
C...     RUNNING SUM
         SUM=SUM+CF(I)*(1.0D0-EXPT)*PSII
C...
C...  CONTINUE SUMMATION
1     CONTINUE
C...
C...  INCLUDE EXPONENTIAL FACTOR
2     PHI1=EXPZ*SUM
C...
C...  CALCULATION IS COMPLETE
      RETURN
      END
      DOUBLE PRECISION FUNCTION F1(Z)
C...
C...  F1 COMPUTES THE TRANSCENDENTAL FUNCTION THAT DEFINES THE
C...  EIGENVALUES.  F1 IS CALL BY FUNCTION ROOT VIA SUBROUTINE
C...  EIGVAL TO FIND THE ZEROS OF THE TRANSCENDENTAL FUNCTION.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  COMPUTE FUNCTION TO BE ZEROED
      F1=1.0D0/DTAN(Z)-Z/PE+PE/(4.0D0*Z)
      RETURN
      END
      DOUBLE PRECISION FUNCTION F2(Z)
C...
C...  FUNCTION F2 COMPUTES THE PRODUCT OF TWO EIGENFUNCTIONS.  F2
C...  IS CALLED BY FUNCTION SIMP VIA SUBROUTINE ORTHO.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  EIGENFUNCTION PSI(I)
      PSII=PE/(2.0D0*ZMU(IN))*DSIN(ZMU(IN)*Z)+DCOS(ZMU(IN)*Z)
C...
C...  EIGENFUNCTION PSI(J)
      PSIJ=PE/(2.0D0*ZMU(JN))*DSIN(ZMU(JN)*Z)+DCOS(ZMU(JN)*Z)
C...
C...  PRODUCT OF THE EIGENFUNCTIONS
      F2=PSII*PSIJ
      RETURN
      END
      DOUBLE PRECISION FUNCTION F3(Z)
C...
C...  FUNCTION F3 COMPUTES THE NUMERATOR INTEGRAND OF THE FOURIER
C...  COEFFICIENT.  F3 IS CALLED BY FUNCTION SIMP VIA SUBROUTINE
C...  COEFF.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  EIGENFUNCTION PSI(I)
      PSII=PE/(2.0D0*ZMU(IN))*DSIN(ZMU(IN)*Z)+DCOS(ZMU(IN)*Z)
C...
C...  PRODUCT OF EXPONENTIAL AND EIGENFUNCTION
      F3=DEXP((-PE/2.0D0)*Z)*PSII
      RETURN
      END
      DOUBLE PRECISION FUNCTION F4(Z)
C...
C...  FUNCTION F4 COMPUTES THE DENOMINATOR INTEGRAND OF THE FOURIER
C...  COEFFICIENT.  F4 IS CALLED BY FUNCTION SIMP VIA SUBROUTINE
C...  COEFF.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=81,NA=200)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/   C(NZ)
     +      /F/  CT(NZ)
     +      /S/  CZ(NZ), CZZ(NZ)
     +      /C/      ZL,      ZU,     ZLU,       PE,
     +          ZMU(NA),  CF(NA)
     +      /I/      NS,      NQ,      NP,       IP,
     +            NCASE,      IN,      JN
C...
C...  EIGENFUNCTION PSI(I)
      PSII=PE/(2.0D0*ZMU(IN))*DSIN(ZMU(IN)*Z)+DCOS(ZMU(IN)*Z)
C...
C...  EIGENFUNCTION SQUARED
      F4=PSII**2
      RETURN
      END
      DOUBLE PRECISION FUNCTION ROOT(F,ZN,ZP,EPS)
C...
C...  FUNCTION ROOT COMPUTES THE ROOT ON AN EQUATION (ZERO OF A
C...  FUNCTION) BY THE SECANT METHOD.  THE CODE WAS ORIGINALLY
C...  PROVIDED BY YOUNG-SOON CHOE, AND MODIFIED SLIGHTLY BY W. E.
C...  SCHIESSER.
C...
C...  ARGUMENTS
C...
C...     F      EXTERNAL TO COMPUTE THE FUNCTION TO BE ZEROED (INPUT)
C...
C...     ZN     LOWER BOUND ON THE ROOT (INPUT)
C...
C...     ZP     UPPER BOUND ON THE ROOT (INPUT)
C...
C...     EPS    ABSOLUTE ERROR TOLERANCE FOR THE ROOT (INPUT)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  EXTERNAL FOR THE FUNCTION TO BE ZEROED
      EXTERNAL F
C...
C...  BEGIN SEARCH FOR THE ROOT BY THE SECANT METHOD
      IP=0
      IN=0
C...
C...  INITIAL ESTIMATE OF THE ROOT
2     Z=(ZP+ZN)/2.0D0
C...
C...  COMPUTE FUNCTION TO BE ZEROED
      FZ=F(Z)
C...
C...  CONTINUE SEARCH
      IF(FZ.LE.0.0D0)THEN
         ZN=Z
         FM=FZ
         IN=1
      ELSE
         ZP=Z
         FP=FZ
         IP=1
      END IF
      IF((IP*IN).EQ.0)    GO TO 2
      IF((FP-FM).GT.0.1D0)GO TO 2
C...
C...  UPDATE Z AND COMPUTE FUNCTION
3     SL=(FP-FM)/(ZP-ZN)
      Z=ZP-FP/SL
      FZ=F(Z)
C...
C...  TEST FOR CONVERGENCE
      IF((DABS(FZ).LT.EPS).OR.(DABS(ZP-ZN).LT.EPS))THEN
C...
C...     CONVERGENCE ACHIEVED SO STORE ROOT AND RETURN TO CALLING
C...     PROGRAM
         ROOT=Z
C...
C...     CONTINUE SEARCH
         ELSE
         IF(FZ.GE.0.0D0)THEN
            FM=FZ
            ZN=Z
            GO TO 3
         ELSE
            FP=FZ
            ZP=Z
            GO TO 3
            END IF
         END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION SIMP(F,A,B,NQ)
C...
C...  FUNCTION SIMP PERFORMS A NUMERICAL QUADRATURE (INTEGRATION) BY
C...  SIMPSON'S RULE
C...
C...  ARGUMENTS
C...
C...     F      EXTERNAL TO COMPUTE THE INTEGRAND (INPUT)
C...
C...     A      LOWER LIMIT OF THE INTEGRAL (INPUT)
C...
C...     B      UPPER LIMIT OF THE INTEGRAL (INPUT)
C...
C...     NQ     NUMBER OF PANELS USED IN SIMPSON'S RULE (MUST BE
C...            EVEN) (INPUT)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  EXTERNAL FUNCTION FOR INTEGRAND
      EXTERNAL F
C...
C...  CHECK FOR AN EVEN NUMBER OF PANELS
      IF((NQ/2*2).NE.NQ)THEN
         WRITE(*,2)NQ
2        FORMAT(//,' NQ = ',I3,' (NOT EVEN)')
         SIMP=0.0D0
         RETURN
      END IF
C...
C...  INTEGRATION INTERVAL
      DZ=(B-A)/DFLOAT(NQ)
C...
C...  INITIALIZE SUM
      SIMP=(F(A)+4.0D0*F(A+DZ)+F(B))
C...
C...  CONTINUE SUM
      IF(NQ.GT.2)THEN
         DO 1 I=2,NQ-1,2
            SIMP=SIMP
     +          +2.0D0*F(A+DFLOAT(I  )*DZ)
     +          +4.0D0*F(A+DFLOAT(I+1)*DZ)
1        CONTINUE
      END IF
C...
C...  INTEGRAL
      SIMP=(DZ/3.0D0)*SIMP
      RETURN
      END
