!   D02TYF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d02tyfe_mod

!     D02TYF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: mmax = 2, neq = 1, nin = 5,      &
                                              nlbc = 1, nmesh_out = 11,        &
                                              nout = 6, nrbc = 1
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: a
!     .. Local Arrays ..
      Integer                              :: m(1) = (/2/)
    Contains
      Subroutine ffun(x,y,neq,m,f)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: f(neq)
        Real (Kind=nag_wp), Intent (In)      :: y(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: sqrt
!       .. Executable Statements ..
        If (y(1,0)<=0.0E0_nag_wp) Then
          f(1) = 0.0_nag_wp
        Else
          f(1) = (y(1,0))**1.5_nag_wp/sqrt(x)
        End If
        Return
      End Subroutine ffun
      Subroutine fjac(x,y,neq,m,dfdy)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: dfdy(neq,neq,0:*)
        Real (Kind=nag_wp), Intent (In)      :: y(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: sqrt
!       .. Executable Statements ..
        If (y(1,0)<=0.0E0_nag_wp) Then
          dfdy(1,1,0) = 0.0_nag_wp
        Else
          dfdy(1,1,0) = 1.5_nag_wp*sqrt(y(1,0))/sqrt(x)
        End If
        Return
      End Subroutine fjac
      Subroutine gafun(ya,neq,m,nlbc,ga)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: ga(nlbc)
        Real (Kind=nag_wp), Intent (In)      :: ya(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        ga(1) = ya(1,0) - 1.0_nag_wp
        Return
      End Subroutine gafun
      Subroutine gbfun(yb,neq,m,nrbc,gb)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: gb(nrbc)
        Real (Kind=nag_wp), Intent (In)      :: yb(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        gb(1) = yb(1,0)
        Return
      End Subroutine gbfun
      Subroutine gajac(ya,neq,m,nlbc,dgady)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: dgady(nlbc,neq,0:*)
        Real (Kind=nag_wp), Intent (In)      :: ya(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        dgady(1,1,0) = 1.0_nag_wp
        Return
      End Subroutine gajac
      Subroutine gbjac(yb,neq,m,nrbc,dgbdy)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: dgbdy(nrbc,neq,0:*)
        Real (Kind=nag_wp), Intent (In)      :: yb(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        dgbdy(1,1,0) = 1.0_nag_wp
        Return
      End Subroutine gbjac
      Subroutine guess(x,neq,m,y,dym)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: dym(neq)
        Real (Kind=nag_wp), Intent (Inout)   :: y(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        y(1,0) = 1.0_nag_wp - x/a
        y(1,1) = -1.0_nag_wp/a
        dym(1) = 0.0_nag_wp
        Return
      End Subroutine guess
    End Module d02tyfe_mod
    Program d02tyfe

!     D02TYF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02tkf, d02tvf, d02tyf, d02tzf, nag_wp
      Use d02tyfe_mod, Only: a, ffun, fjac, gafun, gajac, gbfun, gbjac, guess, &
                             m, mmax, neq, nin, nlbc, nmesh_out, nout, nrbc
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: ainc, ermx, x
      Integer                              :: i, iermx, ifail, ijermx, liwork, &
                                              lrwork, mxmesh, ncol, nmesh
      Logical                              :: failed
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: mesh(:), rwork(:), tol(:), y(:,:)
      Integer, Allocatable                 :: ipmesh(:), iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02TYF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) ncol, nmesh, mxmesh
      liwork = mxmesh*(11*neq+6)
      lrwork = mxmesh*(109*neq**2+78*neq+7)
      Allocate (mesh(mxmesh),tol(neq),rwork(lrwork),y(neq,0:mmax-1), &
        ipmesh(mxmesh),iwork(liwork))

      Read (nin,*) a
      Read (nin,*) tol(1:neq)

      ainc = a/real(nmesh-1,kind=nag_wp)
      mesh(1) = 0.0E0_nag_wp
      Do i = 2, nmesh - 1
        mesh(i) = mesh(i-1) + ainc
      End Do
      mesh(nmesh) = a

      ipmesh(1) = 1
      ipmesh(2:nmesh-1) = 2
      ipmesh(nmesh) = 1

!     Initialize
      ifail = 0
      Call d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rwork, &
        lrwork,iwork,liwork,ifail)

      Write (nout,99999) tol(1), a

!     Solve
      ifail = -1
      Call d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rwork,iwork,ifail)

      failed = ifail /= 0

!     Extract mesh.
      ifail = -1
      Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rwork,iwork, &
        ifail)

      If (ifail/=1) Then
!       Print mesh statistics
        Write (nout,99998) nmesh, ermx, iermx, ijermx
        Write (nout,99997)(i,ipmesh(i),mesh(i),i=1,nmesh)
      End If
      If (.Not. failed) Then
!       Print solution on output mesh.
        Write (nout,99996)
        x = 0.0_nag_wp
        ainc = a/real(nmesh_out-1,kind=nag_wp)
        Do i = 1, nmesh_out
          ifail = 0
          Call d02tyf(x,y,neq,mmax,rwork,iwork,ifail)
          Write (nout,99995) x, y(1,0:1)
          x = x + ainc
        End Do
      End If

99999 Format (//' Tolerance = ',E8.1,' A = ',F8.2)
99998 Format (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, &
        '  in interval ',I4,' for component ',I4/)
99997 Format (/' Mesh points:'/4(I4,'(',I1,')',E11.4))
99996 Format (/' Computed solution'/'       x     solution   derivative')
99995 Format (' ',F8.2,2F11.5)

    End Program d02tyfe
