!   D02BJF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d02bjfe_mod

!      Data for D02BJF example program

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: n = 3, nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: h, xend
       INTEGER, SAVE                       :: k
!   n: number of differential equations
    CONTAINS
!dec$ attributes offload : mic :: output
       SUBROUTINE output(xsol,y)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: xsol
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: y(*)
!         .. Local Scalars ..
          INTEGER                             :: j
!         .. Intrinsic Functions ..
          INTRINSIC                              real
!         .. Executable Statements ..
          WRITE (nout,99999) xsol, (y(j),j=1,n)
          xsol = xend - real(k,kind=nag_wp)*h
          k = k - 1
          RETURN

99999     FORMAT (1X,F8.2,3F13.4)
!dec$ attributes offload : mic :: output
       END SUBROUTINE output
!dec$ attributes offload : mic :: fcn
       SUBROUTINE fcn(x,y,f)

!         .. 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
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: fcn
       END SUBROUTINE fcn
!dec$ attributes offload : mic :: g
       FUNCTION g(x,y)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: g
       END FUNCTION g
    END MODULE d02bjfe_mod
    PROGRAM d02bjfe

!      D02BJF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02bjf, d02bjw, d02bjx, nag_wp
       USE d02bjfe_mod, ONLY : fcn, g, h, k, n, nin, nout, output, xend
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: tol, x, xinit
       INTEGER                             :: i, icase, ifail, iw, j, kinit
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: w(:), y(:), yinit(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D02BJF 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
       DO icase = 1, 4
          WRITE (nout,*)
          SELECT CASE (icase)
          CASE (1)
             WRITE (nout,99995) icase, 'intermediate output, root-finding'
          CASE (2)
             WRITE (nout,99995) icase, 'no intermediate output, root-finding'
          CASE (3)
             WRITE (nout,99995) icase, 'intermediate output, no root-finding'
          CASE (4)
             WRITE (nout,99995) icase, 'no intermediate output, &
                &no root-finding ( integrate to XEND)'
          END SELECT
          DO j = 4, 5
             tol = 10.0E0_nag_wp**(-j)
             WRITE (nout,*)
             WRITE (nout,99999) ' Calculation with TOL =', tol
             x = xinit
             y(1:n) = yinit(1:n)
             IF (icase/=2) THEN
                WRITE (nout,*) '     X         Y(1)         Y(2)         Y(3)'
                k = kinit
                h = (xend-x)/real(k+1,kind=nag_wp)
             END IF
!            ifail: behaviour on error exit   
!                   =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
             ifail = 0
             SELECT CASE (icase)
             CASE (1)
                CALL d02bjf(x,xend,n,y,fcn,tol,'Default',output,g,w,ifail)
                WRITE (nout,99998) '  Root of Y(1) = 0.0 at', x
                WRITE (nout,99997) '  Solution is', (y(i),i=1,n)
             CASE (2)
                CALL d02bjf(x,xend,n,y,fcn,tol,'Default',d02bjx,g,w,ifail)
                WRITE (nout,99998) '  Root of Y(1) = 0.0 at', x
                WRITE (nout,99997) '  Solution is', (y(i),i=1,n)
             CASE (3)
                CALL d02bjf(x,xend,n,y,fcn,tol,'Default',output,d02bjw,w, &
                   ifail)
             CASE (4)
                WRITE (nout,99996) x, (y(i),i=1,n)
                CALL d02bjf(x,xend,n,y,fcn,tol,'Default',d02bjx,d02bjw,w, &
                   ifail)
                WRITE (nout,99996) x, (y(i),i=1,n)
             END SELECT
          END DO
          IF (icase<4) THEN
             WRITE (nout,*)
          END IF
       END DO

99999  FORMAT (1X,A,E8.1)
99998  FORMAT (1X,A,F7.3)
99997  FORMAT (1X,A,3F13.4)
99996  FORMAT (1X,F8.2,3F13.4)
99995  FORMAT (1X,'Case ',I1,': ',A)
    END PROGRAM d02bjfe
