```!   D05BEF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.
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
!     .. Accessibility Statements ..
Private
Public                           :: cf1, cf2, cg1, cg2, ck1, ck2, sol1,  &
sol2
!     .. 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 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 d05befe_mod, Only: cf1, cf2, cg1, cg2, ck1, ck2, iorder, lct, lwk,   &
nmesh, nout, sol1, sol2
Use nag_library, Only: d05bef, nag_wp, x02ajf
!     .. 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
```