C      ALGORITHM 688, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 17, NO. 2, PP. 153-166.  JUNE, 1991.

C
C     THIS DOUBLE PRECISION DRIVER MAY BE USED TO RUN THE SECOND EXAMPLE
C     IN SECTION 7 OF SINCOVEC AND MADSEN, TOMS, 5, 1979.
C     THE INPUT REQUIRED CONSISTS OF :
C                 NINT  :  THE NUMBER OF SUBINTERVALS TO BE USED IN THE
C                          VARIABLE.
C                 KORD  :  THE ORDER OF THE B-SPLINES TO BE USED ( = DEG
C     SEVERAL SETS OF INPUT VALUES MAY BE GIVEN IN FORMAT 2I4.
C
      DOUBLE PRECISION XLEFT,DTUSED,XBKPT(121),U(1,121),SCTCH(121),
     *                 WORK(8000),T0,TOUT,DT,EPS,DX,PI,EXACTV,ERROR
      INTEGER NQ,NSTEPS,NF,NJ,IWORK(1000),NPDE,NINT,NPTS,KORD,NCC,
     *        MF,I,K,INDEX
      COMMON /ENDPT/ XLEFT
      COMMON /GEAR0/ DTUSED,NQ,NSTEPS,NF,NJ
      COMMON /PIVAL/ PI
      PI = ACOS(-1.0D+0)
      NPDE = 1
      NCC = 2
    1 CONTINUE
      write(*,*) 'hi'
      READ(5,9100,END=9000)NINT,KORD
      NPTS = NINT + 1
      T0 = 0.0
      TOUT = 1.0D+0
      DT = 0.1E-8
      EPS = 1.E-4
      MF = 21
      INDEX = 1
      IWORK(1) = 8000
      IWORK(2) = 1000
      DX = 1.0/DBLE(NPTS-1)
      DO 10 I = 1,NPTS
        XBKPT(I) = DBLE(I-1) * DX
   10 CONTINUE
      XLEFT = XBKPT(1)
      WRITE(6,9200) NINT, KORD, EPS
   20 CONTINUE
      CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX,
     *            WORK,IWORK)
      IF ( INDEX .NE. 0 ) GO TO 70
      WRITE(6,9150) TOUT,DTUSED,NSTEPS
      CALL VALUES(XBKPT,U,SCTCH,NPDE,NPTS,NPTS,1,WORK)
      ERROR = 0.D0
      DO 60 K = 1,NPDE
        DO 41 I = 1,NPTS
          ERROR = ERROR + (U(K,I)-EXACTV(XBKPT(I),TOUT))**2
   41   CONTINUE
   60 CONTINUE
      WRITE(6,9300)SQRT(ERROR)
      TOUT = TOUT * 10.0
      IF ( TOUT .LT. 1.0) GO TO 20
   70 CONTINUE
      WRITE(6,9400) INDEX
      GO TO 1
 9000 STOP
 9100 FORMAT(2I4)
 9150 FORMAT(' T = ',E10.3,' DT = ',E10.3,' TOTAL STEPS = ',I6)
 9200 FORMAT(' NINT = ',I3, '    KORD = ', I3, '    EPS = ', E10.3)
 9300 FORMAT(' L2 ERROR = ', D24.15)
 9400 FORMAT('  INDEX = ',I3)
      END
      SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE)
C
C     THIS IS THE USER SUPPLIED SUBROUTINE TO SPECIFY THE DIFFERENTIAL
C     EQUATIONS.
C
      DOUBLE PRECISION U(NPDE),UX(NPDE),UXX(NPDE),FVAL(NPDE),T,X
      INTEGER NPDE
      DOUBLE PRECISION PI
      COMMON/PIVAL/PI
      FVAL(1) = UXX(1) + PI * PI * DSIN( PI * X )
      RETURN
      END
      SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE)
C
C     THIS ROUTINE SPECIFIES THE BOUNDARY CONDITION EQUATIONS.
C
      INTEGER NPDE
      DOUBLE PRECISION T,X,U(NPDE),UX(NPDE),DZDT(NPDE),
     *                 DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE),
     *                 XLEFT
      COMMON /ENDPT/ XLEFT
      IF ( X .NE. XLEFT ) GO TO 10
        DBDU(1,1) = 1.0
        DBDUX(1,1) = 0.0
        DZDT(1) = 0.0
      RETURN
   10 CONTINUE
        DBDU(1,1) = 1.0
        DBDUX(1,1) = 0.0
        DZDT(1) = 0.0
      RETURN
      END
      SUBROUTINE UINIT(X,U,NPDE)
C
C     UINIT GIVES THE INITIAL CONDITIONS AT T=T0.
C
      INTEGER NPDE
      DOUBLE PRECISION X,U(NPDE)
      U(1) = 1.0
      RETURN
      END
      SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE)
C
C     THIS IS THE OPTIONAL ROUTINE PROVIDED IF THE USER WISHES TO
C     SUPPLY AN ANALYTIC JACOBIAN.
C
      INTEGER NPDE
      DOUBLE PRECISION T,X,U(NPDE),UX(NPDE),UXX(NPDE),
     *       DFDU(NPDE,NPDE),DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE)
      DFDU(1,1) = 0.0
      DFDUX(1,1) = 0.0
      DFDUXX(1,1) = 1.0
      RETURN
      END
      DOUBLE PRECISION FUNCTION EXACTV(X,T)
C
C     THIS FUNCTION PROVIDES THE EXACT ANSWER
C     TO THE DIFFERENTIAL EQUATION.
C
      DOUBLE PRECISION X,T
      DOUBLE PRECISION PI
      COMMON/PIVAL/PI
      EXACTV = 1.0 + SIN(PI*X)*(1.0 - EXP(-PI*PI*T))
      RETURN
      END