```!   D02NGF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

Module d02ngfe_mod

!     D02NGF 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, neq = 3,       &
nin = 5, nout = 6
Integer, Parameter, Public       :: nrw = 50 + 4*neq
Integer, Parameter, Public       :: nwkjac = neq*(neq+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)
r(2) = -ydot(2)
r(3) = -ydot(3)

If (ires==1) Then
r(1) = -alpha*y(1) + beta*y(2)*y(3) + 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 d02ngfe_mod

Program d02ngfe

!     D02NGF Example Main Program

!     .. Use Statements ..
Use nag_library, Only: d02nby, d02ngf, d02ngz, d02nsf, d02nvf, d02nyf,   &
d02xjf, nag_wp, x04abf
Use d02ngfe_mod, Only: iset, itrace, ldysav, neq, nin, nout, nrw,        &
nwkjac, resid
!     .. Implicit None Statement ..
Implicit None
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: h, h0, hmax, hmin, hu, t, tcrit,     &
tcur, tolsf, tout, xout
Integer                          :: i, ifail, imxer, iout, 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(:), sol(:),  &
wkjac(:), y(:), ydot(:), ysav(:,:)
Real (Kind=nag_wp)               :: con(6)
Integer                          :: inform(23)
Logical, Allocatable             :: algequ(:)
Logical                          :: lderiv(2)
!     .. Intrinsic Procedures ..
Intrinsic                        :: real
!     .. Executable Statements ..
Write (nout,*) 'D02NGF Example Program Results'
!     Skip heading in data file
sdysav = maxord + 1
Allocate (atol(neq),rtol(neq),rwork(nrw),sol(neq),wkjac(nwkjac),y(neq),  &
ydot(neq),ysav(ldysav,sdysav),algequ(neq))
outchn = nout
Call x04abf(iset,outchn)

!     Integrate to tout by overshooting tout in one step mode (itask=2)
!     using B.D.F formulae with a functional iteration method.
!     Default values for the array con are used. Employ vector
!     tolerances and the Jacobian is evaluated internally, if necessary.
!     monitr subroutine replaced by NAG dummy routine D02NBY.
!     Interpolation outside D02NGF using D02XJF.

Read (nin,*) hmin, hmax, h0, tcrit
con(1:6) = 0.0_nag_wp

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02nvf(neq,sdysav,maxord,'Functional-iteration',petzld,con,tcrit,   &
hmin,hmax,h0,maxstp,mxhnil,'Average-L2',rwork,ifail)

!     Linear algebra setup required (in case functional iteration
!     encounters any difficulty).
ifail = 0
Call d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail)

xout = 0.02E0_nag_wp
iout = 1

Write (nout,99993)(i,i=1,neq)
Write (nout,99999) t, (y(i),i=1,neq)

!      Soft fail and error messages only

loop1: Do
ifail = -1
Call d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,      &
ifail)

If (ifail==0) Then

ifail = 0
Call d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter,  &
imxer,algequ,inform,ifail)

loop2:    Do
If (tcur-hu<xout .And. xout<=tcur) Then

!             C0 interpolation
ifail = 0
Call d02xjf(xout,sol,neq,ysav,ldysav,sdysav,neq,tcur,nqu,hu,h,   &
ifail)

Write (nout,99999) xout, (sol(i),i=1,neq)
iout = iout + 1
xout = real(iout,kind=nag_wp)*0.02E0_nag_wp
If (iout>=6) Then
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
Exit loop1
End If
Else
Exit loop2
End If
End Do loop2
Else
Write (nout,*)
Write (nout,99998) 'Exit D02NGF with IFAIL = ', ifail, '  and T = ', &
t
End If
End Do loop1

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,'    X ',3('         Y(',I1,')  '))
End Program d02ngfe
```