!   D02EJF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    Module d02ejfe_mod

!     Data for D02EJF example program

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fcn, g, output, pederv
!     .. 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
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: n = 3, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp), Public, Save :: h, xend
    Contains
      Subroutine fcn(x,y,f)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: f(*)
        Real (Kind=nag_wp), Intent (In) :: y(*)
!       .. Executable Statements ..
        f(1) = -alpha*y(1) + beta*y(2)*y(3)
        f(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2)
        f(3) = gamma*y(2)*y(2)
        Return
      End Subroutine fcn
      Subroutine pederv(x,y,pw)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: pw(*)
        Real (Kind=nag_wp), Intent (In) :: y(*)
!       .. Executable Statements ..
        pw(1) = -alpha
        pw(2) = alpha
        pw(3) = zero
        pw(4) = beta*y(3)
        pw(5) = -beta*y(3) - 2.0_nag_wp*gamma*y(2)
        pw(6) = 2.0_nag_wp*gamma*y(2)
        pw(7) = beta*y(2)
        pw(8) = -beta*y(2)
        pw(9) = zero
        Return
      End Subroutine pederv
      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) - 0.9E0_nag_wp
        Return
      End Function g
      Subroutine output(xsol,y)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: xsol
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: y(*)
!       .. Local Scalars ..
        Integer                        :: j
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs
!       .. Executable Statements ..
        Write (nout,99999) xsol, (y(j),j=1,n)
        xsol = xsol + h
!       Make sure we exactly hit xsol = xend
        If (abs(xsol-xend)<h/4.0E0_nag_wp) Then
          xsol = xend
        End If
        Return

99999   Format (1X,F8.2,3F13.5)
      End Subroutine output
    End Module d02ejfe_mod
    Program d02ejfe

!     D02EJF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02ejf, d02ejw, d02ejx, d02ejy, nag_wp
      Use d02ejfe_mod, Only: fcn, g, h, n, nin, nout, output, pederv, 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 Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02EJF Example Program Results'
      iw = (12+n)*n + 50
      Allocate (w(iw),y(n),yinit(n))
!     Skip heading in data file
      Read (nin,*)
!     xinit: initial x value, xend: final x value
!     y: initial solution values
      Read (nin,*) xinit, xend
      Read (nin,*) yinit(1:n)
      Read (nin,*) kinit
      Do icase = 1, 5
        If (icase/=2) Then
          Write (nout,99995) icase, 'Jacobian internally'
        Else
          Write (nout,99995) icase, 'Jacobian by PEDERV'
        End If
        Select Case (icase)
        Case (1,2)
          Write (nout,99994) 'intermediate output, root-finding'
        Case (3)
          Write (nout,99994) 'no intermediate output, root-finding'
        Case (4)
          Write (nout,99994) 'intermediate output, no root-finding'
        Case (5)
          Write (nout,99994)                                                   &
            'no intermediate output, no root-finding (integrate to XEND)'
        End Select
        Do j = 3, 4
          tol = 10.0E0_nag_wp**(-j)
          Write (nout,99999) ' Calculation with TOL =', tol
          x = xinit
          y(1:n) = yinit(1:n)
          If (icase/=3) Then
            Write (nout,*) '     X         Y(1)         Y(2)         Y(3)'
            h = (xend-x)/real(kinit+1,kind=nag_wp)
          End If
          ifail = 0
          Select Case (icase)
          Case (1)
            Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',output,g,w,iw,     &
              ifail)
            Write (nout,99998) '  Root of Y(1)-0.9 at', x
            Write (nout,99997) '  Solution is', (y(i),i=1,n)
          Case (2)
            Call d02ejf(x,xend,n,y,fcn,pederv,tol,'Default',output,g,w,iw,     &
              ifail)
            Write (nout,99998) '  Root of Y(1)-0.9 at', x
            Write (nout,99997) '  Solution is', (y(i),i=1,n)
          Case (3)
            Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',d02ejx,g,w,iw,     &
              ifail)
            Write (nout,99998) '  Root of Y(1)-0.9 at', x
            Write (nout,99997) '  Solution is', (y(i),i=1,n)
          Case (4)
            ifail = 0
            Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',output,d02ejw,w,   &
              iw,ifail)
          Case (5)
            Write (nout,99996) x, (y(i),i=1,n)
            Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',d02ejx,d02ejw,w,   &
              iw,ifail)
            Write (nout,99996) x, (y(i),i=1,n)
          End Select
          If (tol<0.0E0_nag_wp) Then
            Write (nout,*) '  Range too short for TOL'
          End If
        End Do
        If (icase<5) Then
          Write (nout,*)
        End If
      End Do

99999 Format (/,1X,A,E8.1)
99998 Format (1X,A,F7.3)
99997 Format (1X,A,3F13.5)
99996 Format (1X,F8.2,3F13.5)
99995 Format (/,1X,'Case ',I1,': calculating ',A,',')
99994 Format (8X,A)
    End Program d02ejfe