!   D05BEF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.
    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
!dec$ attributes offload : mic :: sol1
       FUNCTION sol1(t)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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

!dec$ attributes offload : mic :: sol1
       END FUNCTION sol1
!dec$ attributes offload : mic :: sol2
       FUNCTION sol2(t)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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

!dec$ attributes offload : mic :: sol2
       END FUNCTION sol2
!dec$ attributes offload : mic :: ck1
       FUNCTION ck1(t)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: ck1
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          ck1 = exp(-0.5_nag_wp*t)

          RETURN

!dec$ attributes offload : mic :: ck1
       END FUNCTION ck1
!dec$ attributes offload : mic :: cf1
       FUNCTION cf1(t)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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

!dec$ attributes offload : mic :: cf1
       END FUNCTION cf1
!dec$ attributes offload : mic :: cg1
       FUNCTION cg1(s,y)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: cg1
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: s, y
!         .. Executable Statements ..
          cg1 = y

          RETURN

!dec$ attributes offload : mic :: cg1
       END FUNCTION cg1
!dec$ attributes offload : mic :: ck2
       FUNCTION ck2(t)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          INTRINSIC                              sqrt
!         .. Executable Statements ..
          ck2 = sqrt(x01aaf(pi))

          RETURN

!dec$ attributes offload : mic :: ck2
       END FUNCTION ck2
!dec$ attributes offload : mic :: cf2
       FUNCTION cf2(t)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          INTRINSIC                              log, sqrt
!         .. Executable Statements ..
          st1 = sqrt(1.0_nag_wp+t)
          cf2 = -2.0_nag_wp*log(st1+sqrt(t))/st1

          RETURN

!dec$ attributes offload : mic :: cf2
       END FUNCTION cf2
!dec$ attributes offload : mic :: cg2
       FUNCTION cg2(s,y)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: cg2
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: s, y
!         .. Executable Statements ..
          cg2 = y

          RETURN

!dec$ attributes offload : mic :: cg2
       END FUNCTION cg2
    END MODULE d05befe_mod
    PROGRAM d05befe

!      D05BEF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d05bef, x02ajf
       USE d05befe_mod, ONLY : cf1, cf2, cg1, cg2, ck1, ck2, iorder, lct, lwk, &
                               nag_wp, 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 Functions ..
       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
