*     NAG D02MWF Example Program Text
*     Mark 22 Release. NAG Copyright 2008.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          NEQ, MAXORD, LICOM, LCOM
      PARAMETER        (NEQ=5,MAXORD=5,LICOM=50+NEQ,LCOM=40+(MAXORD+4)
     +                 *NEQ+NEQ**2)
      INTEGER          IJAC
      PARAMETER        (IJAC=1)
*     .. Local Scalars ..
      DOUBLE PRECISION G1, G2, H0, HMAX, T, TOUT
      INTEGER          I, IFAIL, ITASK, ITOL, NADV
      CHARACTER*8      JCEVAL
*     .. Local Arrays ..
      DOUBLE PRECISION ATOL(NEQ), COM(LCOM), RTOL(NEQ), RUSER(1),
     +                 Y(NEQ), YDOT(NEQ)
      INTEGER          ICOM(LICOM), IUSER(1)
*     .. External Subroutines ..
      EXTERNAL         D02MWF, D02NEF, JAC, RES, X04ABF
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D02MWF Example Program Results'
      WRITE (NOUT,*)
      NADV = NOUT
      CALL X04ABF(1,NADV)
      T = 0.0D0
      TOUT = 1.0D0
      ITOL = 1
      DO I = 1, NEQ
         RTOL(I) = 1.0D-8
         ATOL(I) = 1.0D-8
      END DO
*     Set initial values
      Y(1) = 1.0D0
      Y(2) = 0.0D0
      Y(3) = 0.0D0
      Y(4) = 1.0D0
      Y(5) = 1.0D0
      IF (IJAC.EQ.1) THEN
         JCEVAL = 'Analytic'
      ELSE
         JCEVAL = 'Numeric'
      END IF
      HMAX = 0.0D0
      H0 = 0.0D0
      IFAIL = 1
      CALL D02MWF(NEQ,MAXORD,JCEVAL,HMAX,H0,ITOL,ICOM,LICOM,COM,LCOM,
     +            IFAIL)
      IF (IFAIL.EQ.0) THEN
         WRITE (NOUT,*)
     +     't           y(1)       y(2)       y(3)       y(4)      y(5)'
         WRITE (NOUT,99998) T, (Y(I),I=1,NEQ)
         DO I = 1, NEQ
            YDOT(I) = 0.0D0
         END DO
         ITASK = 0
   20    CALL D02NEF(NEQ,T,TOUT,Y,YDOT,RTOL,ATOL,ITASK,RES,JAC,ICOM,COM,
     +               LCOM,IUSER,RUSER,IFAIL)
         WRITE (*,99998) T, (Y(I),I=1,NEQ)
         WRITE (NOUT,*)
         WRITE (NOUT,99999) ITASK
         WRITE (NOUT,*)
         IF ((ITASK.GE.0) .AND. (ITASK.LE.3)) THEN
            IF (T.LT.TOUT) GO TO 20
            G1 = Y(1)**2 + Y(2)**2 - 1
            G2 = Y(1)*Y(3) + Y(2)*Y(4)
            WRITE (NOUT,99997) G1
            WRITE (NOUT,99996) G2
         END IF
      ELSE
         WRITE (NOUT,99995) IFAIL
      END IF
*
99999 FORMAT (1X,'D02NEF returned with ITASK = ',I4)
99998 FORMAT (1X,F6.4,5(F11.6))
99997 FORMAT (1X,'The position-level constraint G1 = ',E12.4)
99996 FORMAT (1X,'The velocity-level constraint G2 = ',E12.4)
99995 FORMAT (1X,' ** D02MWF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE RES(NEQ,T,Y,YDOT,R,IRES,IUSER,RUSER)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          IRES, NEQ
*     .. Array Arguments ..
      DOUBLE PRECISION R(NEQ), RUSER(1), Y(NEQ), YDOT(NEQ)
      INTEGER          IUSER(1)
*     .. Executable Statements ..
      R(1) = Y(3) - YDOT(1)
      R(2) = Y(4) - YDOT(2)
      R(3) = -Y(5)*Y(1) - YDOT(3)
      R(4) = -Y(5)*Y(2) - 1.D0 - YDOT(4)
      R(5) = Y(3)**2 + Y(4)**2 - Y(5) - Y(2)
      RETURN
      END
*
      SUBROUTINE JAC(NEQ,T,Y,YDOT,PD,CJ,IUSER,RUSER)
*     .. Scalar Arguments ..
      DOUBLE PRECISION CJ, T
      INTEGER          NEQ
*     .. Array Arguments ..
      DOUBLE PRECISION PD(NEQ,NEQ), RUSER(1), Y(NEQ), YDOT(NEQ)
      INTEGER          IUSER(1)
*     .. Executable Statements ..
      PD(1,1) = -CJ
      PD(1,3) = 1.0D0
      PD(2,2) = -CJ
      PD(2,4) = 1.0D0
      PD(3,1) = -Y(5)
      PD(3,3) = -CJ
      PD(3,5) = -Y(1)
      PD(4,2) = -Y(5)
      PD(4,4) = -CJ
      PD(4,5) = -Y(2)
      PD(5,2) = -1.0D0
      PD(5,3) = 2.0D0*Y(3)
      PD(5,4) = 2.0D0*Y(4)
      PD(5,5) = -1.0D0
      RETURN
      END