*       D02NEF Example Program Text
*       Mark 22 Release. NAG Copyright 2008
*
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. External Subroutines ..
      EXTERNAL         EX1, EX2
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D02NEF Example Program Results'
      CALL EX1
      CALL EX2
      END
*
      SUBROUTINE EX1
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          NEQ, MAXORD, LICOM, MU, ML, LCOM
      PARAMETER        (NEQ=3,MAXORD=5,LICOM=50+NEQ,MU=2,ML=1,
     +                 LCOM=40+(MAXORD+4)*NEQ+(2*ML+MU+1)
     +                 *NEQ+2*(NEQ/(ML+MU+1)+1))
      INTEGER          IJAC
      PARAMETER        (IJAC=1)
*     .. Local Scalars ..
      DOUBLE PRECISION H0, HMAX, T, TOUT
      INTEGER          I, IFAIL, ITASK, ITOL, J
      CHARACTER*8      JCEVAL
*     .. Local Arrays ..
      DOUBLE PRECISION ATOL(NEQ), COM(LCOM), RTOL(NEQ), RUSER(1),
     +                 Y(NEQ), YDOT(NEQ)
      INTEGER          ICOM(LICOM), IUSER(3)
*     .. External Subroutines ..
      EXTERNAL         D02MCF, D02MWF, D02NEF, D02NPF, JAC1, RES1
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'D02NEF Example 1'
      WRITE (NOUT,*)
      ITOL = 1
      DO I = 1, NEQ
         RTOL(I) = 1.0D-3
         ATOL(I) = 1.0D-6
         YDOT(I) = 0.0D0
      END DO
      IF (IJAC.EQ.1) THEN
         JCEVAL = 'Analytic'
      ELSE
         JCEVAL = 'Numeric'
      END IF
*     Set initial values
      Y(1) = 1.0D0
      Y(2) = 0.0D0
      Y(3) = 0.0D0
*
*     Initialize the problem, specifying that the Jacobian is to be
*     evaluated analytically using the provided routine JAC.
*
      HMAX = 0.0D0
      H0 = 0.0D0
      T = 0.0D0
      TOUT = 0.02D0
*
      IFAIL = 1
      CALL D02MWF(NEQ,MAXORD,JCEVAL,HMAX,H0,ITOL,ICOM,LICOM,COM,LCOM,
     +            IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,99996) IFAIL
         GO TO 40
      END IF
*
*     Specify that the Jacobian is banded.
*
      IFAIL = 0
      CALL D02NPF(NEQ,ML,MU,ICOM,LICOM,IFAIL)
*
*     Use the IUSER array to pass the band dimensions through to JAC.
*     An alternative would be to hard code values for ML and MU in JAC.
*
      IUSER(1) = ML
      IUSER(2) = MU
      IUSER(3) = IJAC
*
      WRITE (NOUT,*) '    t            y(1)        y(2)        y(3)   '
      WRITE (NOUT,99998) T, (Y(I),I=1,NEQ)
      ITASK = 0
*
*     Obtain the solution at 5 equally spaced values of T.
*
      DO 20 J = 1, 5
         IFAIL = 1
         CALL D02NEF(NEQ,T,TOUT,Y,YDOT,RTOL,ATOL,ITASK,RES1,JAC1,ICOM,
     +               COM,LCOM,IUSER,RUSER,IFAIL)
         WRITE (NOUT,99998) T, (Y(I),I=1,NEQ)
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99997) IFAIL
            GO TO 40
         END IF
         TOUT = TOUT + 0.02D0
         CALL D02MCF(ICOM)
   20 CONTINUE
*
      WRITE (NOUT,*)
      WRITE (NOUT,99999) ITASK
*
   40 CONTINUE
*
99999 FORMAT (1X,'The integrator completed task, ITASK = ',I4)
99998 FORMAT (1X,F8.4,3X,3(F12.6))
99997 FORMAT (1X,' ** D02NEF returned with IFAIL = ',I5)
99996 FORMAT (1X,' ** D02MWF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE RES1(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) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3) - YDOT(1)
      R(2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2)*Y(2) - YDOT(2)
      R(3) = 3.0D7*Y(2)*Y(2) - YDOT(3)
      RETURN
      END
*
      SUBROUTINE JAC1(NEQ,T,Y,YDOT,PD,CJ,IUSER,RUSER)
*     .. Scalar Arguments ..
      DOUBLE PRECISION CJ, T
      INTEGER          NEQ
*     .. Array Arguments ..
      DOUBLE PRECISION PD(*), RUSER(1), Y(NEQ), YDOT(NEQ)
      INTEGER          IUSER(3)
*     .. Local Scalars ..
      INTEGER          IJAC, ML, MU
*     .. External Subroutines ..
      EXTERNAL         D02NEZ, MYJAC1
*     .. Executable Statements ..
      ML = IUSER(1)
      MU = IUSER(2)
      IJAC = IUSER(3)
*
      IF (IJAC.EQ.1) THEN
         CALL MYJAC1(NEQ,ML,MU,T,Y,YDOT,PD,CJ)
      ELSE
         CALL D02NEZ(NEQ,T,Y,YDOT,PD,CJ,IUSER,RUSER)
      END IF
*
      RETURN
      END
*
      SUBROUTINE MYJAC1(NEQ,ML,MU,T,Y,YDOT,PD,CJ)
*     .. Scalar Arguments ..
      DOUBLE PRECISION CJ, T
      INTEGER          ML, MU, NEQ
*     .. Array Arguments ..
      DOUBLE PRECISION PD(2*ML+MU+1,NEQ), Y(NEQ), YDOT(NEQ)
*     .. Local Scalars ..
      INTEGER          MD, MS
*     .. Executable Statements ..
*     Main diagonal PDFULL(i,i), i=1,NEQ
      MD = MU + ML + 1
      PD(MD,1) = -0.04D0 - CJ
      PD(MD,2) = -1.0D4*Y(3) - 6.0D7*Y(2) - CJ
      PD(MD,3) = -CJ
*     1 Sub-diagonal PDFULL(i+1:i), i=1,NEQ-1
      MS = MD + 1
      PD(MS,1) = 0.04D0
      PD(MS,2) = 6.0D7*Y(2)
*     First super-diagonal PDFULL(i-1,i), i=2, NEQ
      MS = MD - 1
      PD(MS,2) = 1.0D4*Y(3)
      PD(MS,3) = -1.0D4*Y(2)
*     Second super-diagonal PDFULL(i-2,i), i=3, NEQ
      MS = MD - 2
      PD(MS,3) = 1.0D4*Y(2)
*
      RETURN
      END
*
      SUBROUTINE EX2
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          NEQ, MAXORD, LICOM, LCOM
      PARAMETER        (NEQ=1,MAXORD=5,LICOM=50+NEQ,LCOM=40+(MAXORD+4)
     +                 *NEQ+NEQ*NEQ)
      INTEGER          IJAC
      PARAMETER        (IJAC=1)
*     .. Local Scalars ..
      DOUBLE PRECISION H0, HMAX, T, TOUT
      INTEGER          I, IFAIL, ITASK, ITOL, J
      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         D02MCF, D02MWF, D02NEF, JAC2, RES2
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'D02NEF Example 2'
      WRITE (NOUT,*)
      ITOL = 1
      DO I = 1, NEQ
         RTOL(I) = 0.0D0
         ATOL(I) = 1.0D-8
         YDOT(I) = 0.0D0
      END DO
      IF (IJAC.EQ.1) THEN
         JCEVAL = 'Analytic'
      ELSE
         JCEVAL = 'Numeric'
      END IF
*
*     Initialize the problem, specifying that the Jacobian is to be
*     evaluated analytically using the provided routine JAC.
*
      HMAX = 0.0D0
      H0 = 0.0D0
      T = 0.0D0
      TOUT = 0.2D0
      Y(1) = 2.0D0
*
      IFAIL = 1
      CALL D02MWF(NEQ,MAXORD,JCEVAL,HMAX,H0,ITOL,ICOM,LICOM,COM,LCOM,
     +            IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,99996) IFAIL
         GO TO 40
      END IF
*
*     Use the IUSER array to pass whether numerical or analytic Jacobian
*     is to be used.
*
      IUSER(1) = IJAC
*
      WRITE (NOUT,*) '    t            y(1)'
      WRITE (NOUT,99998) T, (Y(I),I=1,NEQ)
      ITASK = 0
*
*     Obtain the solution at 5 equally spaced values of T.
*
      DO 20 J = 1, 5
         IFAIL = 1
         CALL D02NEF(NEQ,T,TOUT,Y,YDOT,RTOL,ATOL,ITASK,RES2,JAC2,ICOM,
     +               COM,LCOM,IUSER,RUSER,IFAIL)
         WRITE (NOUT,99998) T, (Y(I),I=1,NEQ)
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99997) IFAIL
            GO TO 40
         END IF
         TOUT = TOUT + 0.2D0
         CALL D02MCF(ICOM)
   20 CONTINUE
*
      WRITE (NOUT,*)
      WRITE (NOUT,99999) ITASK
*
   40 CONTINUE
*
99999 FORMAT (1X,'The integrator completed task, ITASK = ',I4)
99998 FORMAT (1X,F8.4,3X,3(F12.6))
99997 FORMAT (1X,' ** D02NEF returned with IFAIL = ',I5)
99996 FORMAT (1X,' ** D02MWF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE RES2(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)
*     .. Intrinsic Functions ..
      INTRINSIC        EXP
*     .. Executable Statements ..
      R(1) = 4.0D0 - Y(1)**2 + T*0.1D0*EXP(Y(1))
      RETURN
      END
*
      SUBROUTINE JAC2(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)
*     .. Local Scalars ..
      INTEGER          IJAC
*     .. External Subroutines ..
      EXTERNAL         D02NEZ, MYJAC2
*     .. Executable Statements ..
      IJAC = IUSER(1)
*
      IF (IJAC.EQ.1) THEN
         CALL MYJAC2(NEQ,T,Y,YDOT,PD,CJ)
      ELSE
         CALL D02NEZ(NEQ,T,Y,YDOT,PD,CJ,IUSER,RUSER)
      END IF
*
      RETURN
      END
*
      SUBROUTINE MYJAC2(NEQ,T,Y,YDOT,PD,CJ)
*     .. Scalar Arguments ..
      DOUBLE PRECISION CJ, T
      INTEGER          NEQ
*     .. Array Arguments ..
      DOUBLE PRECISION PD(NEQ*NEQ), Y(NEQ), YDOT(NEQ)
*     .. Intrinsic Functions ..
      INTRINSIC        EXP
*     .. Executable Statements ..
      PD(1) = -2.0D0*Y(1) + 0.1D0*T*Y(1)*EXP(Y(1))
*
      RETURN
      END