!   D02BJ_P0W_F Example Program Text
!   Mark 28.4 Release. NAG Copyright 2022.

    Module d02bj_p0w_fe_mod

!     Data for D02BJ_P0W_F example program

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      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(ad_handle,x,y,f,iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: f(*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: y(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: alpha, beta
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, tan
!       .. Executable Statements ..
        alpha = ruser(1)
        beta = ruser(2)
        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
      Subroutine g(ad_handle,x,y,retval,iuser,ruser)

!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Real (Kind=nag_wp), Intent (Out) :: retval
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: y(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..
        retval = y(1)
        Return
      End Subroutine g
    End Module d02bj_p0w_fe_mod
    Program d02bj_p0w_fe

!     D02BJ_P0W_F Example Main Program

!     .. Use Statements ..
      Use d02bj_p0w_fe_mod, Only: fcn, g, n, nin, nout
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: d02bj_p0w_f, d02bj_p0w_x
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: alpha = -0.032E0_nag_wp
      Real (Kind=nag_wp), Parameter    :: beta = -0.02E0_nag_wp
!     .. Local Scalars ..
      Type (c_ptr)                     :: ad_handle
      Real (Kind=nag_wp)               :: tol, x, xend, xinit
      Integer                          :: i, ifail, iw, kinit
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: ruser(4)
      Real (Kind=nag_wp), Allocatable  :: w(:), y(:), yinit(:)
      Integer                          :: iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02BJ_P0W_F Example Program Results'

      iw = 20*n
      Allocate (w(iw),y(n),yinit(n))
!     Skip heading in data file
      Read (nin,*)
!     xinit: initial x value,        xend: final x value.
!     yinit: initial solution values
      Read (nin,*) xinit, xend
      Read (nin,*) yinit(1:n)
      Read (nin,*) kinit

      Write (nout,99996) 'no intermediate output, root-finding'
      tol = 1.0E-5_nag_wp
      Write (nout,*)
      Write (nout,99999) ' Calculation with TOL =', tol

      ruser(1) = alpha
      ruser(2) = beta
      ruser(3) = (xend-xinit)/real(kinit+1,kind=nag_wp)
      ruser(4) = xend

      x = xinit
      y(1:n) = yinit(1:n)
      ifail = 0
      Call d02bj_p0w_f(ad_handle,x,xend,n,y,fcn,tol,'Default',d02bj_p0w_x,g,w, &
        iuser,ruser,ifail)

      Write (nout,99998) '  Root of Y(1) = 0.0 at', x
      Write (nout,99997) '  Solution is', (y(i),i=1,n)

99999 Format (1X,A,E8.1)
99998 Format (1X,A,F7.3)
99997 Format (1X,A,3F13.4)
99996 Format (1X,'Case : ',A)

    End Program d02bj_p0w_fe