!   D02BHF Example Program Text
!   Mark 29.3 Release. NAG Copyright 2023.

    Module d02bhfe_mod

!     D02BHF 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                           :: fcn, g
!     .. Parameters ..
      Integer, Parameter, Public       :: n = 3, nin = 5, nout = 6
!     n: number of differential equations
    Contains
      Subroutine fcn(x,y,f)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: alpha = -0.032E0_nag_wp
        Real (Kind=nag_wp), Parameter  :: beta = -0.02E0_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(*)
        Real (Kind=nag_wp), Intent (In) :: y(*)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, tan
!       .. Executable Statements ..
        f(1) = tan(y(3))
        f(2) = alpha*tan(y(3))/y(2) + beta*y(2)/cos(y(3))
        f(3) = alpha/y(2)**2
        Return
      End Subroutine fcn

      Function g(x,y)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: g
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: y(*)
!       .. Executable Statements ..
        g = y(1)
        Return
      End Function g
    End Module d02bhfe_mod

    Program d02bhfe

!     D02BHF Example Main Program

!     .. Use Statements ..
      Use d02bhfe_mod, Only: fcn, g, n, nin, nout
      Use nag_library, Only: d02bhf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: hmax, tol, x, xend, xinit
      Integer                          :: i, ifail, irelab, j
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: w(:,:), y(:), yinit(:)
!     .. Executable Statements ..
      Write (nout,*) 'D02BHF Example Program Results'
      Allocate (w(n,7),y(n),yinit(n))
!     Skip heading in data file
      Read (nin,*)
!     xinit: initial x value,         xend  : final x value.
!     yinit: initial solution values, irelab: type of error control.
      Read (nin,*) xinit
      Read (nin,*) xend
      Read (nin,*) yinit(1:n)
      Read (nin,*) irelab
      hmax = 0.0E0_nag_wp
      Do i = 4, 5
        tol = 10.0E0_nag_wp**(-i)
        x = xinit
        y(1:n) = yinit(1:n)

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d02bhf(x,xend,n,y,tol,irelab,hmax,fcn,g,w,ifail)

        Write (nout,*)
        Write (nout,99999) 'Calculation with TOL =', tol
        Write (nout,99998) ' Root of Y(1) at', x
        Write (nout,99997) ' Solution is', (y(j),j=1,n)
        If (tol<0.0E0_nag_wp) Then
          Write (nout,*) ' Over one-third steps controlled by HMAX'
        End If
      End Do

99999 Format (1X,A,E8.1)
99998 Format (1X,A,F7.4)
99997 Format (1X,A,3F13.5)
    End Program d02bhfe