C
C     THIS DOUBLE PRECISION DRIVER MAY BE USED TO RUN THE FIRST 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(2,121),SCTCH(121),
     *                 WORK(28000),T0,TOUT,DT,EPS,DX
*
      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
      NPDE = 2
    1 CONTINUE
      READ(5,9100,END=9000) NINT, KORD
      NPTS = NINT + 1
      NCC = 2
      T0 = 0.0
      TOUT = 1.D-1
      DT = 1.D-7
      EPS = 1.0D-4
      MF = 21
      WRITE(6,9200)NINT,KORD,EPS
      INDEX = 1
      IWORK(1) = 28000
      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)
   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,0,WORK)
      DO 60 K = 1,NPDE
        WRITE(6,9300)K
        WRITE(6,9400)(U(K,I),I=1,NPTS)
   60 CONTINUE
C     TOUT = TOUT * 10.0
C     IF ( TOUT .LT. 11.0 ) GO TO 20
   70 WRITE(6,9500) INDEX
      GO TO 1
 9000 STOP
 9100 FORMAT(2I4)
 9150 FORMAT(' T = ',E10.3,' DT = ',E10.3,' TOTAL STEPS = ',I6)
 9200 FORMAT('NO. OF SUBINTERVALS = ',I3,'  KORD = ',I2,
     *' EPS = ',D10.2)
 9300 FORMAT(10X,'PDE COMPONENT = ',I3)
 9400 FORMAT(10X,5E12.4)
 9500 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
      FVAL(1) = U(2)*U(2)*UXX(1) - U(1)*U(2) - U(1)*U(1) + 10.0
     *          + 2.0*U(2)*UX(2)*UX(1)
      FVAL(2) = U(1)*U(1)*UXX(2) + U(1)*U(2) - U(2)*U(2)
     *          + UXX(1) + 2.0*U(1)*UX(1)*UX(2)
      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
        DBDU(1,2) = 0.0
        DBDU(2,1) = 0.0
        DBDU(2,2) = 1.0
        DBDUX(1,1) = 0.0
        DBDUX(1,2) = 0.0
        DBDUX(2,1) = 0.0
        DBDUX(2,2) = 0.0
        DZDT(1) = 0.0
        DZDT(2) = 0.0
      RETURN
   10 CONTINUE
        DBDU(1,1) = U(2)*COS(U(1)*U(2))
        DBDU(1,2) = U(1)*COS(U(1)*U(2))
        DBDU(2,1) = U(2)*SIN(U(1)*U(2))
        DBDU(2,2) = U(1)*SIN(U(1)*U(2))
        DBDUX(1,1) = 1.0
        DBDUX(1,2) = 0.0
        DBDUX(2,1) = 0.0
        DBDUX(2,2) = 1.0
        DZDT(1) = 0.0
        DZDT(2) = 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) = 0.5*(X + 1.0)
      U(2) = 4.0*ATAN(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) = -U(2) -2.0*U(1)
      DFDU(1,2) = 2.0*U(2)*UXX(1) - U(1) + 2.0*UX(2)*UX(1)
      DFDU(2,1) = 2.0*U(1)*UXX(2) + U(2) + 2.0*UX(1)*UX(2)
      DFDU(2,2) = U(1) - 2.0*U(2)
      DFDUX(1,1) = 2.0*U(2)*UX(2)
      DFDUX(1,2) = 2.0*U(2)*UX(1)
      DFDUX(2,1) = 2.0*U(1)*UX(2)
      DFDUX(2,2) = 2.0*U(1)*UX(1)
      DFDUXX(1,1) = U(2)*U(2)
      DFDUXX(1,2) = 0.0
      DFDUXX(2,1) = 1.0
      DFDUXX(2,2) = U(1)*U(1)
      RETURN
      END