C An example of the use of DDRIV2
C
C Notice that NROOT is set to 1
C
      DOUBLE PRECISION  H
      PARAMETER (N=2, H=10.0, NROOT=1, MINT=2, 
     *           LW=N*N+10*N+2*NROOT+204, LIW=23)
      DOUBLE PRECISION  Y(N+1),W(LW),GFUN,EPS,T,TOUT,EWT,MASS
      INTEGER   IW(LIW)
      EXTERNAL  FSUB,GFUN
      DATA      MASS /0.125D0/
C
      EPS = 1.D-5
C Set initial point
      T = 0.0D0
      TOUT = T
C Set for pure relative error
      EWT = 0.0D0
C Set initial conditions
      Y(1) = H
      Y(2) = 0.0D0
C Set parameter value
      Y(3) = MASS
      MS = 1
      WRITE (*,*) 'DDRIV2 RESULTS'
      WRITE (*,5) T,Y(1),Y(2),MS
    5 FORMAT(2X,D18.10,2X,D18.10,2X,D18.10,1X,I3)
C
   10 CALL DDRIV2
     *   (N,T,Y,FSUB,TOUT,MS,NROOT,EPS,EWT,MINT,W,LW,IW,LIW,GFUN)
      TOUT = TOUT+0.1D0
      IF (MS .EQ. 5) THEN
        WRITE (*,'(2X,D18.10,2X,D18.10,2X,D18.10,I4)') 
     *        T,Y(1),Y(2),MS
        WRITE (*,'(2X,A,D18.10)') ' <-- Y=0 AT T= ',T
        GOTO 100
      ELSE
        WRITE (*,'(2X,D18.10,2X,D18.10,2X,D18.10,I4)') T,Y(1),Y(2),MS
C Stop if any output code but 1 or 2.
        IF (MS .GT. 2) GOTO 100
      END IF
      GOTO 10
  100 WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS (LAST LINE) FROM IBM PC/AT'
      WRITE (*,*)
     *'   0.2624999995E+01   -0.5648151218E-15   -0.3999999828E+01   5', 
     *'  <-- Y=0 AT T=   0.2624999995E+01'
      STOP
      END 
C
      SUBROUTINE FSUB (N,T,Y,YDOT)
C Routine for evaluating right hand sides of equations.
        DOUBLE PRECISION T, Y(*), YDOT(*), MASS, G
        DATA G/32.0D0/
C Retrieve parameter
        MASS = Y(3)
C
        YDOT(1) = Y(2)
        YDOT(2) = -(G+1.0D0/MASS*Y(2))
        RETURN
      END
C
      DOUBLE PRECISION FUNCTION GFUN(N,T,Y,IROOT)
C Routine for root finding.
C Integration will stop when GFUN changes sign.
        DOUBLE PRECISION Y(*),T
        GFUN = Y(1)
        RETURN
      END