*     D02CJF Example Program Text
*     Mark 14 Revised.  NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          N, IW
      PARAMETER        (N=3,IW=21*N+28)
*     .. 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 D02CJW, G, X01AAF
      EXTERNAL         D02CJW, G, X01AAF
*     .. External Subroutines ..
      EXTERNAL         D02CJF, D02CJX, FCN, OUT
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Common blocks ..
      COMMON           XEND, H, K
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D02CJF 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
         K = 4
         H = (XEND-X)/DBLE(K+1)
         WRITE (NOUT,*) '     X         Y(1)         Y(2)         Y(3)'
         IFAIL = 0
*
         CALL D02CJF(X,XEND,N,Y,FCN,TOL,'Default',OUT,G,W,IFAIL)
*
         WRITE (NOUT,99998) '  Root of Y(1) = 0.0 at', X
         WRITE (NOUT,99997) '  Solution is', (Y(I),I=1,N)
   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 = 0
*
         CALL D02CJF(X,XEND,N,Y,FCN,TOL,'Default',D02CJX,G,W,IFAIL)
*
         WRITE (NOUT,99998) '  Root of Y(1) = 0.0 at', X
         WRITE (NOUT,99997) '  Solution is', (Y(I),I=1,N)
   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 = 0
*
         CALL D02CJF(X,XEND,N,Y,FCN,TOL,'Default',OUT,D02CJW,W,IFAIL)
*
   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 = 0
*
         CALL D02CJF(X,XEND,N,Y,FCN,TOL,'Default',D02CJX,D02CJW,W,IFAIL)
*
         WRITE (NOUT,99996) X, (Y(I),I=1,N)
   80 CONTINUE
      STOP
*
99999 FORMAT (1X,A,D8.1)
99998 FORMAT (1X,A,F7.3)
99997 FORMAT (1X,A,3F13.5)
99996 FORMAT (1X,F8.2,3F13.5)
      END
*
      SUBROUTINE OUT(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.5)
      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