*     D02NMF Example Program Text
*     Mark 14 Revised. NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          NEQ, NEQMAX, NRW, NINF, NWKJAC, NJCPVT, MAXORD,
     +                 SDYSAV, MAXSTP, MXHNIL, LDYSAV
      PARAMETER        (NEQ=3,NEQMAX=NEQ,NRW=50+4*NEQMAX,NINF=23,
     +                 NWKJAC=NEQMAX*(NEQMAX+1),NJCPVT=1,MAXORD=5,
     +                 SDYSAV=MAXORD+1,MAXSTP=200,MXHNIL=5,LDYSAV=NEQ)
      DOUBLE PRECISION H0, HMAX, HMIN, TCRIT
      PARAMETER        (H0=0.0D0,HMAX=10.0D0,HMIN=1.0D-10,TCRIT=0.0D0)
      LOGICAL          PETZLD
      PARAMETER        (PETZLD=.FALSE.)
*     .. Local Scalars ..
      DOUBLE PRECISION H, HLAST, HNEXT, HU, T, TC, TCUR, TOLSF, TOUT,
     +                 XOUT
      INTEGER          I, IFAIL, IFLAG, IMON, IMXER, INLN, IOUT, IRES,
     +                 IREVCM, ITASK, ITOL, ITRACE, LACOR1, LACOR2,
     +                 LACOR3, LACORB, LSAVR1, LSAVR2, LSAVR3, LSAVRB,
     +                 NITER, NJE, NQ, NQU, NRE, NST, OUTCHN
*     .. Local Arrays ..
      DOUBLE PRECISION ATOL(NEQ), CONST(6), RTOL(NEQ), RWORK(NRW),
     +                 WKJAC(NWKJAC), Y(NEQ), YDOT(NEQ),
     +                 YSAV(LDYSAV,SDYSAV)
      INTEGER          INFORM(NINF), JACPVT(NJCPVT)
      LOGICAL          ALGEQU(NEQ)
*     .. External Subroutines ..
      EXTERNAL         D02NMF, D02NSF, D02NVF, D02NYF, D02XKF, X04ABF
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE, INT
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D02NMF Example Program Results'
      OUTCHN = NOUT
      WRITE (NOUT,*)
      CALL X04ABF(1,OUTCHN)
*
*     Integrate to TOUT by overshooting TOUT (ITASK=1) using B.D.F.
*     formulae with a Newton method. Default values for the array CONST
*     are used. Employ scalar tolerances and the Jacobian is evaluated
*     internally. On the reverse communication call equivalent to the
*     MONITR call in forward communication routines carry out
*     interpolation using D02XKF.
*
      T = 0.0D0
      TOUT = 10.0D0
      ITASK = 1
      IOUT = 1
      XOUT = 2.0D0
      Y(1) = 1.0D0
      Y(2) = 0.0D0
      Y(3) = 0.0D0
      ITOL = 1
      RTOL(1) = 1.0D-4
      ATOL(1) = 1.0D-7
      DO 20 I = 1, 6
         CONST(I) = 0.0D0
   20 CONTINUE
      IFAIL = 1
*
      CALL D02NVF(NEQMAX,SDYSAV,MAXORD,'Newton',PETZLD,CONST,TCRIT,HMIN,
     +            HMAX,H0,MAXSTP,MXHNIL,'Average-L2',RWORK,IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,99993) IFAIL
         GO TO 80
      END IF
      IFAIL = 0
      CALL D02NSF(NEQ,NEQMAX,'Numerical',NWKJAC,RWORK,IFAIL)
*
      LACORB = 50 + NEQ
      LACOR1 = LACORB + 1
      LACOR2 = LACORB + 2
      LACOR3 = LACORB + 3
      LSAVRB = LACORB + NEQ
      LSAVR1 = LSAVRB + 1
      LSAVR2 = LSAVRB + 2
      LSAVR3 = LSAVRB + 3
      WRITE (NOUT,*) '    X          Y(1)           Y(2)           Y(3)'
      WRITE (NOUT,99999) T, (Y(I),I=1,NEQ)
*
*     Soft fail and error messages only
      IREVCM = 0
      ITRACE = 0
   40 IFAIL = 1
*
      CALL D02NMF(NEQ,LDYSAV,T,TOUT,Y,YDOT,RWORK,RTOL,ATOL,ITOL,INFORM,
     +            YSAV,SDYSAV,WKJAC,NWKJAC,JACPVT,NJCPVT,IMON,INLN,IRES,
     +            IREVCM,ITASK,ITRACE,IFAIL)
*
      IF (IREVCM.NE.0) THEN
         IF (IREVCM.EQ.1 .OR. IREVCM.EQ.3) THEN
*           Equivalent to FCN evaluation in forward communication
*           routines
            RWORK(LSAVR1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3)
            RWORK(LSAVR2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2)
     +                      *Y(2)
            RWORK(LSAVR3) = 3.0D7*Y(2)*Y(2)
         ELSE IF (IREVCM.EQ.4) THEN
*           Equivalent to FCN evaluation in forward communication
*           routines
            RWORK(LACOR1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3)
            RWORK(LACOR2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2)
     +                      *Y(2)
            RWORK(LACOR3) = 3.0D7*Y(2)*Y(2)
         ELSE IF (IREVCM.EQ.5) THEN
*           Equivalent to FCN evaluation in forward communication
*           routines
            YDOT(1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3)
            YDOT(2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2)*Y(2)
            YDOT(3) = 3.0D7*Y(2)*Y(2)
         ELSE IF (IREVCM.EQ.9) THEN
*           Equivalent to MONITR call in forward communication routines
            IF (IMON.EQ.1) THEN
               TC = RWORK(19)
               HLAST = RWORK(15)
               HNEXT = RWORK(16)
               NQU = INT(RWORK(10))
   60          CONTINUE
               IF (TC-HLAST.LT.XOUT .AND. XOUT.LE.TC) THEN
                  IFLAG = 1
*
                  CALL D02XKF(XOUT,RWORK(LSAVR1),NEQ,YSAV,LDYSAV,SDYSAV,
     +                        RWORK(LACOR1),NEQ,TC,NQU,HLAST,HNEXT,
     +                        IFLAG)
*
                  IF (IFLAG.NE.0) THEN
                     IMON = -2
                  ELSE
                     WRITE (NOUT,99999) XOUT, (RWORK(LSAVRB+I),I=1,NEQ)
                     IOUT = IOUT + 1
                     XOUT = DBLE(IOUT)*2.0D0
                     IF (IOUT.LT.6) GO TO 60
                  END IF
               END IF
            END IF
         ELSE IF (IREVCM.EQ.2 .OR. IREVCM.EQ.6 .OR. IREVCM.EQ.7 .OR.
     +            IREVCM.EQ.8) THEN
            WRITE (NOUT,*)
            WRITE (NOUT,99994) 'Illegal value of IREVCM = ', IREVCM
            GO TO 80
         END IF
         GO TO 40
      ELSE
         IF (IFAIL.EQ.0) THEN
*
            IFAIL = 0
            CALL D02NYF(NEQ,NEQMAX,HU,H,TCUR,TOLSF,RWORK,NST,NRE,NJE,
     +                  NQU,NQ,NITER,IMXER,ALGEQU,INFORM,IFAIL)
*
            WRITE (NOUT,*)
            WRITE (NOUT,99997) HU, H, TCUR
            WRITE (NOUT,99996) NST, NRE, NJE
            WRITE (NOUT,99995) NQU, NQ, NITER
            WRITE (NOUT,99994) ' Max err comp = ', IMXER
            WRITE (NOUT,*)
         ELSE
            WRITE (NOUT,*)
            WRITE (NOUT,99998) 'Exit D02NMF with IFAIL = ', IFAIL,
     +        '  and T = ', T
         END IF
      END IF
   80 CONTINUE
*
99999 FORMAT (1X,F8.3,3(F13.5,2X))
99998 FORMAT (1X,A,I2,A,E12.5)
99997 FORMAT (1X,' HUSED = ',E12.5,'  HNEXT = ',E12.5,'  TCUR = ',E12.5)
99996 FORMAT (1X,' NST = ',I6,'    NRE = ',I6,'    NJE = ',I6)
99995 FORMAT (1X,' NQU = ',I6,'    NQ  = ',I6,'  NITER = ',I6)
99994 FORMAT (1X,A,I4)
99993 FORMAT (1X,/1X,' ** D02NVF returned with IFAIL = ',I5)
      END