*     D02ZAF Example Program Text
*     Mark 22 Release. NAG Copyright 2008.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          NEQ, NEQMAX, NRW, NINF, NWKJAC, MAXORD, SDYSAV,
     +                 MAXSTP, MXHNIL, LDYSAV
      PARAMETER        (NEQ=3,NEQMAX=NEQ,NRW=50+4*NEQMAX,NINF=23,
     +                 NWKJAC=NEQMAX*(NEQMAX+1),MAXORD=5,
     +                 SDYSAV=MAXORD+1,MAXSTP=200,MXHNIL=5,LDYSAV=NEQ)
      DOUBLE PRECISION H0, HMAX, HMIN
      PARAMETER        (H0=0.0D0,HMAX=10.0D0,HMIN=1.0D-10)
      LOGICAL          PETZLD
      PARAMETER        (PETZLD=.FALSE.)
*     .. Local Scalars ..
      DOUBLE PRECISION H, HU, T, TCRIT, TCUR, TOLSF, TOUT
      INTEGER          I, IFAIL, IMXER, ITASK, ITOL, ITRACE, 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)
      LOGICAL          ALGEQU(NEQ)
*     .. External Subroutines ..
      EXTERNAL         D02NBF, D02NSF, D02NVF, D02NYF, FCN, JAC, MONITR,
     +                 X04ABF
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D02ZAF Example Program Results'
      OUTCHN = NOUT
      CALL X04ABF(1,OUTCHN)
*
*     Set problem parameters
*
      T = 0.0D0
      TOUT = 10.0D0
*
*              Initialisation
*     Integrate to TOUT without passing TOUT.
*
      TCRIT = TOUT
      ITASK = 4
*
*     Use default values for the array CONST.
*
      DO 20 I = 1, 6
         CONST(I) = 0.0D0
   20 CONTINUE
*
*     Use BDF formulae with modified Newton method.
*     Use averaged L2 norm for local error control.
*
      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 40
      END IF
*
*     Setup for using analytic Jacobian
*
      IFAIL = 0
      CALL D02NSF(NEQ,NEQMAX,'Analytical',NWKJAC,RWORK,IFAIL)
*
      WRITE (NOUT,*)
      WRITE (NOUT,*) '  Analytic Jacobian'
      WRITE (NOUT,*)
*
*     Set tolerances.
*
      ITOL = 1
      RTOL(1) = 1.0D-4
      ATOL(1) = 1.0D-7
*
*     Initial values for Y.
*
      Y(1) = 1.0D0
      Y(2) = 0.0D0
      Y(3) = 0.0D0
      WRITE (NOUT,*) '    X          Y(1)           Y(2)           Y(3)'
      WRITE (NOUT,99999) T, (Y(I),I=1,NEQ)
*
*     Solve the problem using MONITR to output intermediate results.
*
      ITRACE = 0
      IFAIL = 1
*
      CALL D02NBF(NEQ,LDYSAV,T,TOUT,Y,YDOT,RWORK,RTOL,ATOL,ITOL,INFORM,
     +            FCN,YSAV,SDYSAV,JAC,WKJAC,NWKJAC,MONITR,ITASK,ITRACE,
     +            IFAIL)
*
      IF (IFAIL.EQ.0) THEN
*
*        Get integration statistics.
*
         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 D02NBF with IFAIL = ', IFAIL,
     +     '  and T = ', T
      END IF
   40 CONTINUE
*
99999 FORMAT (1X,F10.6,3(F13.7,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
*
      SUBROUTINE FCN(NEQ,T,Y,R,IRES)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          IRES, NEQ
*     .. Array Arguments ..
      DOUBLE PRECISION R(NEQ), Y(NEQ)
*     .. Executable Statements ..
      R(1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3)
      R(2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2)*Y(2)
      R(3) = 3.0D7*Y(2)*Y(2)
      RETURN
      END
*
      SUBROUTINE JAC(NEQ,T,Y,H,D,P)
*     .. Scalar Arguments ..
      DOUBLE PRECISION D, H, T
      INTEGER          NEQ
*     .. Array Arguments ..
      DOUBLE PRECISION P(NEQ,NEQ), Y(NEQ)
*     .. Local Scalars ..
      DOUBLE PRECISION HXD
*     .. Executable Statements ..
      HXD = H*D
      P(1,1) = 1.0D0 - HXD*(-0.04D0)
      P(1,2) = -HXD*(1.0D4*Y(3))
      P(1,3) = -HXD*(1.0D4*Y(2))
      P(2,1) = -HXD*(0.04D0)
      P(2,2) = 1.0D0 - HXD*(-1.0D4*Y(3)-6.0D7*Y(2))
      P(2,3) = -HXD*(-1.0D4*Y(2))
*     Do not need to set P(3,1) since Jacobian preset to zero
*     P(3,1) =       - HXD*(0.0E0)
      P(3,2) = -HXD*(6.0D7*Y(2))
      P(3,3) = 1.0D0 - HXD*(0.0D0)
      RETURN
      END
*
      SUBROUTINE MONITR(NEQ,LDYSAV,T,HLAST,HNEXT,Y,YDOT,YSAV,R,ACOR,
     +                  IMON,INLN,HMIN,HMAX,NQU)
*
*     Fortran Library Implementation Wrapper routine for D02NBYN.
*     NAG COPYRIGHT 2004.
*
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. Scalar Arguments ..
      DOUBLE PRECISION HLAST, HMAX, HMIN, HNEXT, T
      INTEGER          IMON, INLN, LDYSAV, NEQ, NQU
*     .. Array Arguments ..
      DOUBLE PRECISION ACOR(NEQ,*), R(*), Y(*), YDOT(*), YSAV(LDYSAV,*)
*     .. External Functions ..
      DOUBLE PRECISION D02ZAF
      EXTERNAL         D02ZAF
*     .. Local Scalars ..
      DOUBLE PRECISION ERRLOC
      INTEGER          I, IFAIL
*     .. Executable Statements ..
      CONTINUE
      IF (IMON.EQ.1) THEN
         IFAIL = 1
         ERRLOC = D02ZAF(NEQ,ACOR(1,2),ACOR(1,1),IFAIL)
         IF (IFAIL.NE.0) THEN
            IMON = -2
         ELSE IF (ERRLOC.GT.5.0D0) THEN
            WRITE (NOUT,99999) T, (Y(I),I=1,NEQ), ERRLOC
         ELSE
            WRITE (NOUT,99998) T, (Y(I),I=1,NEQ)
         END IF
      END IF
*
      RETURN
*
99999 FORMAT (1X,F10.6,3(F13.7,2X),/1X,' ** WARNING scaled local error',
     +       ' = ',F13.5)
99998 FORMAT (1X,F10.6,3(F13.7,2X))
      END