!   D05BEF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.
    Module d05befe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: iorder = 4
      Integer, Parameter                   :: nmesh = 2**6 + 2*iorder - 1
      Integer, Parameter                   :: nout = 6
      Integer, Parameter                   :: lct = nmesh/32 + 1
      Integer, Parameter                   ::                                  &
                     lwk = (2*iorder+6)*nmesh + 8*iorder*iorder - 16*iorder + 1
    Contains
      Function sol1(t)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Function Return Value ..
        Real (Kind=nag_wp)                   :: sol1
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: c, pi, t1
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp, sqrt
!       .. Executable Statements ..
        t1 = 1.0_nag_wp + t
        c = 1.0_nag_wp/sqrt(2.0_nag_wp*x01aaf(pi))
        sol1 = c*(1.0_nag_wp/(t**1.5_nag_wp))*exp(-t1*t1/(2.0_nag_wp*t))

        Return

      End Function sol1
      Function sol2(t)

!       .. Function Return Value ..
        Real (Kind=nag_wp)                   :: sol2
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
!       .. Executable Statements ..
        sol2 = 1.0_nag_wp/(1.0_nag_wp+t)

        Return

      End Function sol2
      Function ck1(t)

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

        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)                   :: a, pi, t1
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp, sqrt
!       .. Executable Statements ..
        t1 = 1.0_nag_wp + t
        a = 1.0_nag_wp/sqrt(x01aaf(pi)*t)
        cf1 = -a*exp(-0.5_nag_wp*t1*t1/t)

        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

        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
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: st1
!       .. Intrinsic Procedures ..
        Intrinsic                            :: log, sqrt
!       .. Executable Statements ..
        st1 = sqrt(1.0_nag_wp+t)
        cf2 = -2.0_nag_wp*log(st1+sqrt(t))/st1

        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
!       .. Executable Statements ..
        cg2 = y

        Return

      End Function cg2
    End Module d05befe_mod
    Program d05befe

!     D05BEF Example Main Program

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

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

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

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

      errmax = 0.0_nag_wp

      Do i = 2, nmesh
        hi1 = real(i-1,kind=nag_wp)*h
        err = abs(yn(i)-sol1(hi1))

        If (err>errmax) Then
          errmax = err
          t = hi1
          soln = yn(i)
        End If

        If (i>5 .And. mod(i,5)==1) Then
          Write (nout,99998) hi1, yn(i)
        End If

      End Do

      Write (nout,*)
      Write (nout,99999) errmax, t, soln
      Write (nout,*)

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

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

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

      errmax = 0.0_nag_wp

      Do i = 1, nmesh
        hi1 = real(i-1,kind=nag_wp)*h
        err = abs(yn(i)-sol2(hi1))

        If (err>errmax) Then
          errmax = err
          t = hi1
          soln = yn(i)
        End If

        If (i>7 .And. mod(i,7)==1) Then
          Write (nout,99998) hi1, yn(i)
        End If

      End Do

      Write (nout,*)
      Write (nout,99999) errmax, t, soln

99999 Format (' The maximum absolute error, ',E10.2,', occurred at T =', &
        F8.4/' with solution ',F8.4)
99998 Format (1X,F8.4,F15.4)
99997 Format (' The stepsize h = ',F8.4)
    End Program d05befe