dhstrt Subroutine

   subroutine dhstrt(me, neq, a, b, y, yprime, etol, morder, small, big, spy, pv, yp, sf, h)

      implicit none

      class(ddeabm_class), intent(inout)  :: me
      integer, intent(in)                 :: neq       !! the number of (first order) differential equations to be integrated.
      real(wp), intent(in)                :: a         !! the initial point of integration.
      real(wp), intent(in)                :: b         !! a value of the independent variable used to define
                                                    !! the direction of integration. a reasonable choice is to
                                                    !! set `b` to the first point at which a solution is desired.
                                                    !! you can also use `b`, if necessary, to restrict the length
                                                    !! of the first integration step because the algorithm will
                                                    !! not compute a starting step length which is bigger than
                                                    !! `abs(b-a)`, unless `b` has been chosen too close to `a`.
                                                    !! (it is presumed that dhstrt has been called with `b`
                                                    !! different from `a` on the machine being used. also see the
                                                    !! discussion about the parameter `small`.)
      real(wp), dimension(neq), intent(in) :: y         !! the vector of initial values of the `neq` solution
                                                    !! components at the initial point `a`.
      real(wp), dimension(neq), intent(in) :: yprime    !! the vector of derivatives of the `neq`
                                                    !! solution components at the initial point `a`.
                                                    !! (defined by the differential equations in subroutine `me%df`)
      real(wp), dimension(neq), intent(in) :: etol      !! the vector of error tolerances corresponding to
                                                    !! the `neq` solution components. it is assumed that all
                                                    !! elements are positive. following the first integration
                                                    !! step, the tolerances are expected to be used by the
                                                    !! integrator in an error test which roughly requires that
                                                    !! `abs(local error) <= etol` for each vector component.
      integer, intent(in) :: morder    !! the order of the formula which will be used by
                                    !! the initial value method for taking the first integration
                                    !! step.
      real(wp), intent(in) :: small    !! a small positive machine dependent constant
                                    !! which is used for protecting against computations with
                                    !! numbers which are too small relative to the precision of
                                    !! floating point arithmetic.  `small` should be set to
                                    !! (approximately) the smallest positive double precision
                                    !! number such that `(1.0_wp+small) > 1.0_wp`.  on the machine being
                                    !! used. the quantity `small**(3.0_wp/8.0_wp)` is used in computing
                                    !! increments of variables for approximating derivatives by
                                    !! differences.  also the algorithm will not compute a
                                    !! starting step length which is smaller than
                                    !! `100.0_wp*small*abs(a)`.
      real(wp), intent(in) :: big      !! a large positive machine dependent constant which
                                    !! is used for preventing machine overflows. a reasonable
                                    !! choice is to set big to (approximately) the square root of
                                    !! the largest double precision number which can be held in
                                    !! the machine.
      real(wp), intent(out) :: h       !! appropriate starting step size to be attempted by the
                                    !! differential equation method.
      real(wp), dimension(neq), intent(inout) :: spy  !! work array of length `neq` which provide the routine with needed storage space.
      real(wp), dimension(neq), intent(inout) :: pv   !! work array of length `neq` which provide the routine with needed storage space.
      real(wp), dimension(neq), intent(inout) :: yp   !! work array of length `neq` which provide the routine with needed storage space.
      real(wp), dimension(neq), intent(inout) :: sf   !! work array of length `neq` which provide the routine with needed storage space.

      integer :: j, k, lk
      real(wp) :: absdx, da, delf, dely, &
                  dfdub, dfdxb, &
                  dx, dy, fbnd, relper, &
                  srydpb, tolexp, tolmin, tolp, tolsum, ydpb

      if (me%error) return !user-triggered error

      ! begin block permitting
      dx = b - a
      absdx = abs(dx)
      relper = small**0.375_wp

      ! compute an approximate bound (dfdxb) on the partial
      ! derivative of the equation with respect to the
      ! independent variable. protect against an overflow.
      ! also compute a bound (fbnd) on the first derivative
      ! locally.

      da = sign(max(min(relper*abs(a), absdx), 100.0_wp*small*abs(a)), dx)
      if (da == 0.0_wp) da = relper*dx
      call me%df(a + da, y, sf)
      if (me%error) return ! user-triggered error
      yp = sf - yprime
      delf = dhvnrm(yp, neq)
      dfdxb = big
      if (delf < big*abs(da)) dfdxb = delf/abs(da)
      fbnd = dhvnrm(sf, neq)

      ! compute an estimate (dfdub) of the local lipschitz
      ! constant for the system of differential equations. this
      ! also represents an estimate of the norm of the jacobian
      ! locally.  three iterations (two when neq=1) are used to
      ! estimate the lipschitz constant by numerical differences.
      ! the first perturbation vector is based on the initial
      ! derivatives and direction of integration. the second
      ! perturbation vector is formed using another evaluation of
      ! the differential equation.  the third perturbation vector
      ! is formed using perturbations based only on the initial
      ! values. components that are zero are always changed to
      ! non-zero values (except on the first iteration). when
      ! information is available, care is taken to ensure that
      ! components of the perturbation vector have signs which are
      ! consistent with the slopes of local solution curves.
      ! also choose the largest bound (fbnd) for the first
      ! derivative.
      !
      ! perturbation vector size is held
      ! constant for all iterations. compute
      ! this change from the
      ! size of the vector of initial
      ! values.

      dely = relper*dhvnrm(y, neq)
      if (dely == 0.0_wp) dely = relper
      dely = sign(dely, dx)
      delf = dhvnrm(yprime, neq)
      fbnd = max(fbnd, delf)
      if (delf == 0.0_wp) then
         ! cannot have a null perturbation vector
         spy = 0.0_wp
         yp = 1.0_wp
         delf = dhvnrm(yp, neq)
      else
         ! use initial derivatives for first perturbation
         spy = yprime
         yp = yprime
      end if

      dfdub = 0.0_wp
      lk = min(neq + 1, 3)
      do k = 1, lk
         ! define perturbed vector of initial values
         pv = y + yp*(dely/delf)
         if (k == 2) then
            ! use a shifted value of the independent variable
            ! in computing one estimate
            call me%df(a + da, pv, yp)
            if (me%error) return ! user-triggered error
            pv = yp - sf
         else
            ! evaluate derivatives associated with perturbed
            ! vector and compute corresponding differences
            call me%df(a, pv, yp)
            if (me%error) return ! user-triggered error
            pv = yp - yprime
         end if
         ! choose largest bounds on the first derivative
         ! and a local lipschitz constant
         fbnd = max(fbnd, dhvnrm(yp, neq))
         delf = dhvnrm(pv, neq)
         if (delf >= big*abs(dely)) then
            ! protect against an overflow
            dfdub = big
            exit
         end if
         dfdub = max(dfdub, delf/abs(dely))
         if (k == lk) exit

         ! choose next perturbation vector
         if (delf == 0.0_wp) delf = 1.0_wp
         do j = 1, neq
            if (k == 2) then
               dy = y(j)
               if (dy == 0.0_wp) dy = dely/relper
            else
               dy = abs(pv(j))
               if (dy == 0.0_wp) dy = delf
            end if
            if (spy(j) == 0.0_wp) spy(j) = yp(j)
            if (spy(j) /= 0.0_wp) dy = sign(dy, spy(j))
            yp(j) = dy
         end do
         delf = dhvnrm(yp, neq)
      end do

      ! compute a bound (ydpb) on the norm of the second derivative

      ydpb = dfdxb + dfdub*fbnd

      ! define the tolerance parameter upon which the starting step
      ! size is to be based.  a value in the middle of the error
      ! tolerance range is selected.

      tolmin = big
      tolsum = 0.0_wp
      do k = 1, neq
         tolexp = log10(etol(k))
         tolmin = min(tolmin, tolexp)
         tolsum = tolsum + tolexp
      end do
      tolp = 10.0_wp**(0.5_wp*(tolsum/neq + tolmin)/(morder + 1))

      ! compute a starting step size based on the above first and
      ! second derivative information
      !
      ! restrict the step length to be not bigger
      ! than abs(b-a). (unless b is too close to a)

      h = absdx

      if (ydpb == 0.0_wp .and. fbnd == 0.0_wp) then
         ! both first derivative term (fbnd) and second
         ! derivative term (ydpb) are zero
         if (tolp < 1.0_wp) h = absdx*tolp
      else if (ydpb == 0.0_wp) then
         ! only second derivative term (ydpb) is zero
         if (tolp < fbnd*absdx) h = tolp/fbnd
      else
         ! second derivative term (ydpb) is non-zero
         srydpb = sqrt(0.5_wp*ydpb)
         if (tolp < srydpb*absdx) h = tolp/srydpb
      end if

      ! further restrict the step length to be not bigger than  1/dfdub
      if (h*dfdub > 1.0_wp) h = 1.0_wp/dfdub

      ! finally, restrict the step length to be not
      ! smaller than 100*small*abs(a). however, if
      ! a=0. and the computed h underflowed to zero,
      ! the algorithm returns small*abs(b) for the step length.
      h = max(h, 100.0_wp*small*abs(a))
      if (h == 0.0_wp) h = small*abs(b)

      ! now set direction of integration
      h = sign(h, dx)

   end subroutine dhstrt