dvhin Subroutine

   subroutine dvhin(me,n,t0,y0,ydot,tout,uround,ewt,itol,&
                    atol,y,temp,h0,niter,ier)

      class(dvode_t),intent(inout) :: me
      real(wp),intent(in) :: t0 !! initial value of independent variable
      real(wp),intent(in) :: y0(*) !! vector of initial conditions
      real(wp),intent(in) :: ydot(*) !! vector of initial first derivatives
      real(wp),intent(in) :: tout !! first output value of independent variable
      real(wp),intent(in) :: uround !! machine unit roundoff
      real(wp),intent(in) :: ewt(*) !! error weights and tolerance parameters as described in the driver routine
      real(wp),intent(in) :: atol(*) !! error weights and tolerance parameters as described in the driver routine
      real(wp),intent(inout) :: y(n) !! work array of length n
      real(wp),intent(inout) :: temp(n) !! work array of length n
      integer,intent(in) :: n !! size of ode system
      integer,intent(in) :: itol !! error weights and tolerance parameters as described in the driver routine
      real(wp),intent(out) :: h0 !! step size to be attempted
      integer,intent(out) :: niter !! number of iterations (and of f evaluations) to compute h0
      integer,intent(out) :: ier !! the error flag, returned with the value:
                                 !!
                                 !!  * ier = 0  if no trouble occurred, or
                                 !!  * ier = -1 if tout and t0 are considered too close to proceed.

      real(wp) :: afi , atoli , delyi , h , hg , hlb , &
                  hnew , hrat , hub , t1 , tdist , &
                  tround , yddnrm
      integer :: i , iter

      real(wp),parameter :: half = 0.5_wp
      real(wp),parameter :: hun = 100.0_wp
      real(wp),parameter :: pt1 = 0.1_wp
      real(wp),parameter :: two = 2.0_wp

      main : block

         niter = 0
         tdist = abs(tout-t0)
         tround = uround*max(abs(t0),abs(tout))
         if ( tdist<two*tround ) then
            ! error return for tout - t0 too small. --------------------------------
            ier = -1
            return
         else
            ! set a lower bound on h based on the roundoff level in t0 and tout. ---
            hlb = hun*tround
            ! set an upper bound on h based on tout-t0 and the initial y and ydot. -
            hub = pt1*tdist
            atoli = atol(1)
            do i = 1 , n
               if ( itol==2 .or. itol==4 ) atoli = atol(i)
               delyi = pt1*abs(y0(i)) + atoli
               afi = abs(ydot(i))
               if ( afi*hub>delyi ) hub = delyi/afi
            enddo
            ! set initial guess for h as geometric mean of upper and lower bounds. -
            iter = 0
            hg = sqrt(hlb*hub)
            ! if the bounds have crossed, exit with the mean value. ----------------
            if ( hub<hlb ) then
               h0 = hg
               exit main
            endif
            do
               ! looping point for iteration. -----------------------------------------
               ! estimate the second derivative as a difference quotient in f. --------
               h = sign(hg,tout-t0)
               t1 = t0 + h
               do i = 1 , n
                  y(i) = y0(i) + h*ydot(i)
               enddo
               call me%f(n,t1,y,temp)
               do i = 1 , n
                  temp(i) = (temp(i)-ydot(i))/h
               enddo
               yddnrm = me%dvnorm(n,temp,ewt(1:n))
               ! get the corresponding new value of h. --------------------------------
               if ( yddnrm*hub*hub>two ) then
                  hnew = sqrt(two/yddnrm)
               else
                  hnew = sqrt(hg*hub)
               endif
               iter = iter + 1
               !-----------------------------------------------------------------------
               ! test the stopping conditions.
               ! stop if the new and previous h values differ by a factor of < 2.
               ! stop if four iterations have been done.  also, stop with previous h
               ! if hnew/hg > 2 after first iteration, as this probably means that
               ! the second derivative value is bad because of cancellation error.
               !-----------------------------------------------------------------------
               if ( iter<4 ) then
                  hrat = hnew/hg
                  if ( (hrat<=half) .or. (hrat>=two) ) then
                     if ( (iter>=2) .and. (hnew>two*hg) ) then
                        hnew = hg
                        exit
                     endif
                     hg = hnew
                     cycle
                  endif
               endif
               exit
            end do
            ! iteration done.  apply bounds, bias factor, and sign.  then exit. ----
            h0 = hnew*half
            if ( h0<hlb ) h0 = hlb
            if ( h0>hub ) h0 = hub
         endif

      end block main

      h0 = sign(h0,tout-t0)
      niter = iter
      ier = 0

   end subroutine dvhin