!   D02NHF Example Program Text
!   Mark 28.3 Release. NAG Copyright 2022.

    Module d02nhfe_mod

!     D02NHF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: resid
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: alpha = 0.04_nag_wp
      Real (Kind=nag_wp), Parameter    :: beta = 1.0E4_nag_wp
      Real (Kind=nag_wp), Parameter    :: gamma = 3.0E7_nag_wp
      Integer, Parameter, Public       :: iset = 1, itrace = 0, ml = 1,        &
                                          mu = 2, neq = 3, nin = 5
      Integer, Parameter, Public       :: njcpvt = neq
      Integer, Parameter, Public       :: nout = 6
      Integer, Parameter, Public       :: nrw = 50 + 4*neq
      Integer, Parameter, Public       :: nwkjac = neq*(2*ml+mu+1)
      Integer, Parameter, Public       :: ldysav = neq
    Contains
      Subroutine resid(neq,t,y,ydot,r,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: r(neq)
        Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
!       .. Executable Statements ..
        r(1) = -ydot(1) - ydot(2) - ydot(3)
        r(2) = -ydot(2)
        r(3) = -ydot(3)
        If (ires==1) Then
          r(1) = r(1)
          r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) + r(2)
          r(3) = gamma*y(2)*y(2) + r(3)
        End If
        Return
      End Subroutine resid
    End Module d02nhfe_mod

    Program d02nhfe

!     D02NHF Example Main Program

!     .. Use Statements ..
      Use d02nhfe_mod, Only: iset, itrace, ldysav, ml, mu, neq, nin, njcpvt,   &
                             nout, nrw, nwkjac, resid
      Use nag_library, Only: d02nby, d02nhf, d02nhz, d02ntf, d02nvf, d02nyf,   &
                             nag_wp, x04abf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: h, h0, hmax, hmin, hu, t, tcrit,     &
                                          tcur, tolsf, tout, tout1
      Integer                          :: i, ifail, imxer, itask, itol,        &
                                          maxord, maxstp, mxhnil, niter, nje,  &
                                          nq, nqu, nre, nst, outchn, sdysav
      Logical                          :: petzld
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: atol(:), rtol(:), rwork(:),          &
                                          wkjac(:), y(:), ydot(:), ysav(:,:)
      Real (Kind=nag_wp)               :: con(6)
      Integer                          :: inform(23)
      Integer, Allocatable             :: jacpvt(:)
      Logical, Allocatable             :: algequ(:)
      Logical                          :: lderiv(2)
!     .. Executable Statements ..
      Write (nout,*) 'D02NHF Example Program Results'
      Flush (nout)

!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) maxord, maxstp, mxhnil
      sdysav = maxord + 1
      Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq),ydot(neq), &
        ysav(ldysav,sdysav),jacpvt(njcpvt),algequ(neq))

      outchn = nout
      Call x04abf(iset,outchn)

!     Set itask=6 to provide initial estimates of solution and its
!     time derivative. Default values for the array con are used.
!     Use the B.D.F. formulae with a Newton method.
!     Employ scalar relative tolerance and vector absolute tolerance.
!     The Jacobian is evaluated internally.
!     monitr subroutine replaced by NAG dummy routine D02NBY.

      Read (nin,*) hmin, hmax, h0, tcrit
      Read (nin,*) petzld
      Read (nin,*) t, tout1
      Read (nin,*) y(1:neq)
      Read (nin,*) lderiv(1:2)
      Read (nin,*) itol
      Read (nin,*) rtol(1)
      Read (nin,*) atol(1:neq)
      itask = 6
      tout = tout1
      con(1:6) = 0.0_nag_wp

      ifail = 0
      Call d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0,    &
        maxstp,mxhnil,'Average-L2',rwork,ifail)

      ifail = 0
      Call d02ntf(neq,neq,'Numerical',ml,mu,nwkjac,njcpvt,rwork,ifail)

!     Hard fail on initial and subsequent tasks.
      ifail = 0
      Call d02nhf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,resid,  &
        ysav,sdysav,d02nhz,wkjac,nwkjac,jacpvt,njcpvt,d02nby,lderiv,itask,     &
        itrace,ifail)

      Write (nout,*)
      Write (nout,99999) ' Initial    Y : ', y(1:neq)
      Write (nout,99999) ' Initial YDOT : ', ydot(1:neq)
      Flush (nout)

!     Use these initial estimates and integrate to tout (overshoot and
!     interpolate)
      itask = 1
      tout = tout1

      ifail = 0
      Call d02nhf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,resid,  &
        ysav,sdysav,d02nhz,wkjac,nwkjac,jacpvt,njcpvt,d02nby,lderiv,itask,     &
        itrace,ifail)

      Write (nout,99993)(i,i=1,neq)
      Write (nout,99998) tout, y(1:neq)

      ifail = 0
      Call d02nyf(neq,neq,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

99999 Format (1X,A,3(F11.4,2X))
99998 Format (1X,F8.3,3(F13.5,2X))
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,'    X ',3('         Y(',I1,')  '))
    End Program d02nhfe