*     D02BJF Example Program Text
*     Mark 18 Release. NAG Copyright 1997.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          N, IW
      PARAMETER        (N=3,IW=20*N)
*     .. Scalars in Common ..
      DOUBLE PRECISION H, XEND
      INTEGER          K
*     .. Local Scalars ..
      DOUBLE PRECISION PI, TOL, X
      INTEGER          I, IFAIL, J
*     .. Local Arrays ..
      DOUBLE PRECISION W(IW), Y(N)
*     .. External Functions ..
      DOUBLE PRECISION D02BJW, G, X01AAF
      EXTERNAL         D02BJW, G, X01AAF
*     .. External Subroutines ..
      EXTERNAL         D02BJF, D02BJX, FCN, OUTPUT
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Common blocks ..
      COMMON           XEND, H, K
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D02BJF Example Program Results'
      XEND = 10.0D0
      PI = X01AAF(0.0D0)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Case 1: intermediate output, root-finding'
      DO 20 J = 4, 5
         TOL = 10.0D0**(-J)
         WRITE (NOUT,*)
         WRITE (NOUT,99999) ' Calculation with TOL =', TOL
         X = 0.0D0
         Y(1) = 0.5D0
         Y(2) = 0.5D0
         Y(3) = PI/5.0D0
         WRITE (NOUT,*) '     X         Y(1)         Y(2)         Y(3)'
         K = 4
         H = (XEND-X)/DBLE(K+1)
         IFAIL = 1
*
         CALL D02BJF(X,XEND,N,Y,FCN,TOL,'Default',OUTPUT,G,W,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
            WRITE (NOUT,99998) '  Root of Y(1) = 0.0 at', X
            WRITE (NOUT,99997) '  Solution is', (Y(I),I=1,N)
         ELSE
            WRITE (NOUT,99995) IFAIL
            GO TO 100
         END IF
   20 CONTINUE
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Case 2: no intermediate output, root-finding'
      DO 40 J = 4, 5
         TOL = 10.0D0**(-J)
         WRITE (NOUT,*)
         WRITE (NOUT,99999) ' Calculation with TOL =', TOL
         X = 0.0D0
         Y(1) = 0.5D0
         Y(2) = 0.5D0
         Y(3) = PI/5.0D0
         IFAIL = 1
*
         CALL D02BJF(X,XEND,N,Y,FCN,TOL,'Default',D02BJX,G,W,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
            WRITE (NOUT,99998) '  Root of Y(1) = 0.0 at', X
            WRITE (NOUT,99997) '  Solution is', (Y(I),I=1,N)
         ELSE
            WRITE (NOUT,99995) IFAIL
            GO TO 100
         END IF
   40 CONTINUE
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Case 3: intermediate output, no root-finding'
      DO 60 J = 4, 5
         TOL = 10.0D0**(-J)
         WRITE (NOUT,*)
         WRITE (NOUT,99999) ' Calculation with TOL =', TOL
         X = 0.0D0
         Y(1) = 0.5D0
         Y(2) = 0.5D0
         Y(3) = PI/5.0D0
         K = 4
         H = (XEND-X)/DBLE(K+1)
         WRITE (NOUT,*) '     X         Y(1)         Y(2)         Y(3)'
         IFAIL = 1
*
         CALL D02BJF(X,XEND,N,Y,FCN,TOL,'Default',OUTPUT,D02BJW,W,IFAIL)
*
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99995) IFAIL
            GO TO 100
         END IF
   60 CONTINUE
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*)
     +'Case 4: no intermediate output, no root-finding ( integrate to XE
     +ND)'
      DO 80 J = 4, 5
         TOL = 10.0D0**(-J)
         WRITE (NOUT,*)
         WRITE (NOUT,99999) ' Calculation with TOL =', TOL
         X = 0.0D0
         Y(1) = 0.5D0
         Y(2) = 0.5D0
         Y(3) = PI/5.0D0
         WRITE (NOUT,*) '     X         Y(1)         Y(2)         Y(3)'
         WRITE (NOUT,99996) X, (Y(I),I=1,N)
         IFAIL = 1
*
         CALL D02BJF(X,XEND,N,Y,FCN,TOL,'Default',D02BJX,D02BJW,W,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
            WRITE (NOUT,99996) X, (Y(I),I=1,N)
         ELSE
            WRITE (NOUT,99995) IFAIL
            GO TO 100
         END IF
   80 CONTINUE
  100 CONTINUE
*
99999 FORMAT (1X,A,E8.1)
99998 FORMAT (1X,A,F7.3)
99997 FORMAT (1X,A,3F13.4)
99996 FORMAT (1X,F8.2,3F13.4)
99995 FORMAT (1X,/1X,' ** D02BJF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE OUTPUT(X,Y)
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          N
      PARAMETER        (N=3)
*     .. Scalar Arguments ..
      DOUBLE PRECISION X
*     .. Array Arguments ..
      DOUBLE PRECISION Y(N)
*     .. Scalars in Common ..
      DOUBLE PRECISION H, XEND
      INTEGER          I
*     .. Local Scalars ..
      INTEGER          J
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Common blocks ..
      COMMON           XEND, H, I
*     .. Executable Statements ..
      WRITE (NOUT,99999) X, (Y(J),J=1,N)
      X = XEND - DBLE(I)*H
      I = I - 1
      RETURN
*
99999 FORMAT (1X,F8.2,3F13.4)
      END
*
      SUBROUTINE FCN(T,Y,F)
*     .. Parameters ..
      INTEGER          N
      PARAMETER        (N=3)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
*     .. Array Arguments ..
      DOUBLE PRECISION F(N), Y(N)
*     .. Intrinsic Functions ..
      INTRINSIC        COS, TAN
*     .. Executable Statements ..
      F(1) = TAN(Y(3))
      F(2) = -0.032D0*TAN(Y(3))/Y(2) - 0.02D0*Y(2)/COS(Y(3))
      F(3) = -0.032D0/Y(2)**2
      RETURN
      END
*
      DOUBLE PRECISION FUNCTION G(T,Y)
*     .. Parameters ..
      INTEGER          N
      PARAMETER        (N=3)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
*     .. Array Arguments ..
      DOUBLE PRECISION Y(N)
*     .. Executable Statements ..
      G = Y(1)
      RETURN
      END