*     D05BEF Example Program Text
*     Mark 16 Release. NAG Copyright 1992.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          IORDER, NMESH, LCT, LWK
      PARAMETER        (IORDER=4,NMESH=2**6+2*IORDER-1,LCT=NMESH/32+1,
     +                 LWK=(2*IORDER+6)*NMESH+8*IORDER*IORDER-
     +                 16*IORDER+1)
*     .. Local Scalars ..
      DOUBLE PRECISION ERR, ERRMAX, H, HI1, SOLN, T, TLIM, TOLNL
      INTEGER          I, IFAIL
*     .. Local Arrays ..
      DOUBLE PRECISION WORK(LWK), YN(NMESH)
      INTEGER          NCT(LCT)
*     .. External Functions ..
      DOUBLE PRECISION CF1, CF2, CG1, CG2, CK1, CK2, SOL1, SOL2, X02AJF
      EXTERNAL         CF1, CF2, CG1, CG2, CK1, CK2, SOL1, SOL2, X02AJF
*     .. External Subroutines ..
      EXTERNAL         D05BEF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS, DBLE, MOD, SQRT
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D05BEF Example Program Results'
      WRITE (NOUT,*)
      IFAIL = 0
      TLIM = 7.0D0
      TOLNL = SQRT(X02AJF())
      H = TLIM/(NMESH-1)
*
      YN(1) = 0.0D0
*
      CALL D05BEF(CK1,CF1,CG1,'Initial',IORDER,TLIM,TOLNL,NMESH,YN,WORK,
     +            LWK,NCT,IFAIL)
*
      WRITE (NOUT,*) 'Example 1'
      WRITE (NOUT,*)
      WRITE (NOUT,99997) H
      WRITE (NOUT,*)
      WRITE (NOUT,*) '     T        Approximate'
      WRITE (NOUT,*) '                Solution '
      WRITE (NOUT,*)
*
      ERRMAX = 0.0D0
      DO 20 I = 2, NMESH
         HI1 = DBLE(I-1)*H
         ERR = ABS(YN(I)-SOL1(HI1))
         IF (ERR.GT.ERRMAX) THEN
            ERRMAX = ERR
            T = HI1
            SOLN = YN(I)
         END IF
         IF (I.GT.5 .AND. MOD(I,5).EQ.1) WRITE (NOUT,99998) HI1,
     +       YN(I)
   20 CONTINUE
      WRITE (NOUT,*)
      WRITE (NOUT,99999) ERRMAX, T, SOLN
*
      WRITE (NOUT,*)
*
      TLIM = 5.0D0
      H = TLIM/(NMESH-1)
      YN(1) = 1.0D0
*
      CALL D05BEF(CK2,CF2,CG2,'Subsequent',IORDER,TLIM,TOLNL,NMESH,YN,
     +            WORK,LWK,NCT,IFAIL)
*
      WRITE (NOUT,*) 'Example 2'
      WRITE (NOUT,*)
      WRITE (NOUT,99997) H
      WRITE (NOUT,*)
      WRITE (NOUT,*) '     T        Approximate'
      WRITE (NOUT,*) '                Solution '
      WRITE (NOUT,*)
*
      ERRMAX = 0.0D0
      DO 40 I = 1, NMESH
         HI1 = DBLE(I-1)*H
         ERR = ABS(YN(I)-SOL2(HI1))
         IF (ERR.GT.ERRMAX) THEN
            ERRMAX = ERR
            T = HI1
            SOLN = YN(I)
         END IF
         IF (I.GT.7 .AND. MOD(I,7).EQ.1) WRITE (NOUT,99998) HI1, YN(I)
   40 CONTINUE
      WRITE (NOUT,*)
      WRITE (NOUT,99999) ERRMAX, T, SOLN
*
      STOP
*
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
*
*
      DOUBLE PRECISION FUNCTION CK1(T)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              T
*     .. Intrinsic Functions ..
      INTRINSIC                     EXP
*     .. Executable Statements ..
      CK1 = EXP(-0.5D0*T)
      RETURN
      END
*
*
      DOUBLE PRECISION FUNCTION CF1(T)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              T
*     .. Local Scalars ..
      DOUBLE PRECISION              A, PI, T1
*     .. External Functions ..
      DOUBLE PRECISION              X01AAF
      EXTERNAL                      X01AAF
*     .. Intrinsic Functions ..
      INTRINSIC                     EXP, SQRT
*     .. Executable Statements ..
      T1 = 1.0D0 + T
      A = 1.0D0/SQRT(X01AAF(PI)*T)
      CF1 = -A*EXP(-0.5D0*T1*T1/T)
      RETURN
      END
*
*
      DOUBLE PRECISION FUNCTION CG1(S,Y)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              S, Y
*     .. Executable Statements ..
      CG1 = Y
      RETURN
      END
*
*
      DOUBLE PRECISION FUNCTION SOL1(T)
*     .. Scalar Arguments ..
      DOUBLE PRECISION               T
*     .. Local Scalars ..
      DOUBLE PRECISION               C, PI, T1
*     .. External Functions ..
      DOUBLE PRECISION               X01AAF
      EXTERNAL                       X01AAF
*     .. Intrinsic Functions ..
      INTRINSIC                      EXP, SQRT
*     .. Executable Statements ..
*
      T1 = 1.0D0 + T
      C = 1.0D0/SQRT(2.0D0*X01AAF(PI))
      SOL1 = C*(1.0D0/(T**1.5D0))*EXP(-T1*T1/(2.0D0*T))
      RETURN
      END
*
*
      DOUBLE PRECISION FUNCTION CK2(T)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              T
*     .. Local Scalars ..
      DOUBLE PRECISION              PI
*     .. External Functions ..
      DOUBLE PRECISION              X01AAF
      EXTERNAL                      X01AAF
*     .. Intrinsic Functions ..
      INTRINSIC                     SQRT
*     .. Executable Statements ..
      CK2 = SQRT(X01AAF(PI))
      RETURN
      END
*
*
      DOUBLE PRECISION FUNCTION CF2(T)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              T
*     .. Local Scalars ..
      DOUBLE PRECISION              ST1
*     .. Intrinsic Functions ..
      INTRINSIC                     LOG, SQRT
*     .. Executable Statements ..
      ST1 = SQRT(1.0D0+T)
      CF2 = -2.0D0*LOG(ST1+SQRT(T))/ST1
      RETURN
      END
*
*
      DOUBLE PRECISION FUNCTION CG2(S,Y)
*     .. Scalar Arguments ..
      DOUBLE PRECISION              S, Y
*     .. Executable Statements ..
      CG2 = Y
      RETURN
      END
*
*
      DOUBLE PRECISION FUNCTION SOL2(T)
*     .. Scalar Arguments ..
      DOUBLE PRECISION               T
*     .. Executable Statements ..
      SOL2 = 1.0D0/(1.0D0+T)
      RETURN
      END