!   D05BDF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.
    Module d05bdfe_mod

!     D05BDF 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                           :: cf1, cf2, cg1, cg2, ck1, ck2
!     .. Parameters ..
      Integer, Parameter, Public       :: iorder = 4
      Integer, Parameter, Public       :: nmesh = 2**6 + 2*iorder - 1
      Integer, Parameter, Public       :: nout = 6
      Integer, Parameter, Public       :: lct = nmesh/32 + 1
      Integer, Parameter, Public       :: lwk = (2*iorder+6)*nmesh + 8*iorder* &
                                          iorder - 16*iorder + 1
    Contains
      Function ck1(t)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: ck1
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: pi
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        ck1 = -sqrt(x01aaf(pi))

        Return

      End Function ck1
      Function cf1(t)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: cf1
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: pi
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        cf1 = sqrt(t) + (3.0_nag_wp/8.0_nag_wp)*t*t*x01aaf(pi)

        Return

      End Function cf1
      Function cg1(s,y)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: cg1
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: s, y
!       .. Executable Statements ..
        cg1 = y*y*y

        Return

      End Function cg1
      Function ck2(t)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: ck2
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: pi
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        ck2 = -sqrt(x01aaf(pi))

        Return

      End Function ck2
      Function cf2(t)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: cf2
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        cf2 = (3.0_nag_wp-t)*sqrt(t)

        Return

      End Function cf2
      Function cg2(s,y)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: cg2
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: s, y
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        cg2 = exp(s*(1.0_nag_wp-s)*(1.0_nag_wp-s)-y*y)

        Return

      End Function cg2
    End Module d05bdfe_mod
    Program d05bdfe

!     D05BDF Example Main Program

!     .. Use Statements ..
      Use d05bdfe_mod, Only: cf1, cf2, cg1, cg2, ck1, ck2, iorder, lct, lwk,   &
                             nmesh, nout
      Use nag_library, Only: d05bdf, nag_wp, x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: h, t, tlim, tolnl
      Integer                          :: i, ifail
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: work(lwk), yn(nmesh)
      Integer                          :: nct(lct)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: mod, real, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'D05BDF Example Program Results'

      tlim = 7.0_nag_wp
      tolnl = sqrt(x02ajf())
      h = tlim/real(nmesh-1,kind=nag_wp)

      ifail = 0
      Call d05bdf(ck1,cf1,cg1,'Initial',iorder,tlim,tolnl,nmesh,yn,work,lwk,   &
        nct,ifail)

      Write (nout,*)
      Write (nout,*) 'Example 1'
      Write (nout,*)
      Write (nout,99998) h
      Write (nout,*)
      Write (nout,*) '     T        Approximate'
      Write (nout,*) '                Solution '
      Write (nout,*)

      Do i = 1, nmesh
        t = real(i-1,kind=nag_wp)*h

        If (mod(i,5)==1) Then
          Write (nout,99999) t, yn(i)
        End If

      End Do

      tlim = 5.0_nag_wp
      h = tlim/real(nmesh-1,kind=nag_wp)

      ifail = 0
      Call d05bdf(ck2,cf2,cg2,'Subsequent',iorder,tlim,tolnl,nmesh,yn,work,    &
        lwk,nct,ifail)

      Write (nout,*)
      Write (nout,*) 'Example 2'
      Write (nout,*)
      Write (nout,99998) h
      Write (nout,*)
      Write (nout,*) '     T        Approximate'
      Write (nout,*) '                Solution '
      Write (nout,*)

      Do i = 1, nmesh
        t = real(i-1,kind=nag_wp)*h

        If (mod(i,7)==1) Then
          Write (nout,99999) t, yn(i)
        End If

      End Do

99999 Format (1X,F8.4,F15.4)
99998 Format (' The stepsize h = ',F8.4)
    End Program d05bdfe