ddeabm_initialize Subroutine

   subroutine ddeabm_initialize(me, neq, maxnum, df, rtol, atol, report, initial_step_mode, initial_step_size)

      implicit none

      class(ddeabm_class), intent(inout)   :: me
      integer, intent(in)                  :: neq     !! number of equations to be integrated
      integer, intent(in)                  :: maxnum  !! maximum number of integration steps allowed
      procedure(deriv_func)               :: df      !! derivative function dx/dt
      real(wp), dimension(:), intent(in)    :: rtol    !! relative tolerance for integration (see [[ddeabm]])
      real(wp), dimension(:), intent(in)    :: atol    !! absolution tolerance for integration (see [[ddeabm]])
      procedure(report_func), optional     :: report  !! reporting function
      integer, intent(in), optional         :: initial_step_mode !! how to choose the initial step `h`:
                                                             !!
                                                             !! 1. Use [[dhstrt]].
                                                             !! 2. Use the older (quicker) algorithm.
                                                             !! 3. Use the user-specified value `initial_step_size` (>0).
      real(wp), intent(in), optional        :: initial_step_size !! for `initial_step_mode=3`

      logical :: vector_tols
      integer :: i, j, k

      integer, parameter :: max_order = 12 !!@Note the max order of this code is 12.
                                           !! This value is hard-coded in several places (array dimensions, etc.)
                                           !! Eventually, may want to make this a user-selectable value.
                                           !! For now, it is fixed.

      !initialize the class:
      call me%destroy()

      !number of equations:
      me%neq = neq

      !maximum number of steps:
      me%maxnum = maxnum

      !set the derivative routine pointer:
      me%df => df

      !set the report function (optional):
      if (present(report)) then
         me%report => report
      else
         nullify (me%report)
      end if

      !allocate the work arrays:
      allocate (me%ypout(neq))
      allocate (me%yp(neq))
      allocate (me%yy(neq))
      allocate (me%wt(neq))
      allocate (me%p(neq))
      allocate (me%phi(neq, 16))

      !tolerances:
      ! [for now, we are considering these unchangeable, although they do not have to be]

      vector_tols = size(rtol) == neq .and. size(atol) == neq .and. neq > 1
      me%scalar_tols = .not. vector_tols

      if (me%scalar_tols) then
         allocate (me%rtol(1)); allocate (me%rtol_tmp(1))
         allocate (me%atol(1)); allocate (me%atol_tmp(1))
         me%rtol = rtol(1)
         me%atol = atol(1)
      else
         allocate (me%rtol(size(rtol))); allocate (me%rtol_tmp(size(rtol)))
         allocate (me%atol(size(atol))); allocate (me%atol_tmp(size(atol)))
         me%rtol = rtol
         me%atol = atol
      end if

      !Compute the two and gstr arrays:
      !  Note: this is a modification of the original dsteps code.
      !  The full-precision coefficients are used here, instead
      !  of the less precise ones in the original.
      !  These are computed recursively using the equation on p. 159 of
      !  Shampine/Gordon, "Computer Solution of Ordinary Differential Equations", 1975.
      allocate (me%two(max_order + 1))
      allocate (me%gstr(0:max_order + 1))
      do i = 1, max_order + 1
         me%two(i) = 2.0_wp**i
      end do
      me%gstr = 0.0_wp
      me%gstr(0) = 1.0_wp
      do i = 1, max_order + 1
         k = 1
         do j = i - 1, 0, -1
            k = k + 1
            me%gstr(i) = me%gstr(i) + 1.0_wp/real(k, wp)*me%gstr(j)
         end do
         me%gstr(i) = -me%gstr(i)
      end do
      me%gstr = abs(me%gstr)

      ! initial step size (can have initial_step_mode 1, 2, or 3):
      me%initial_step_size = 0.0_wp ! dummy value (only used if initial_step_mode is 3)
      if (present(initial_step_mode)) then
         select case (initial_step_mode)
         case (1, 2, 3)
            me%initial_step_mode = initial_step_mode
            if (initial_step_mode == 3) then
               if (present(initial_step_size)) then
                  if (initial_step_size > 0.0_wp) then
                     me%initial_step_size = initial_step_size
                  else
                     error stop 'Error in ddeabm_initialize: initial_step_size must be > 0'
                  end if
               else
                  error stop 'Error in ddeabm_initialize: initial_step_size must be input when initial_step_mode=3'
               end if
            end if
         case default
            error stop 'Error in ddeabm_initialize: invalid input for initial_step_mode'
         end select
      else
         me%initial_step_mode = 1 ! original method (dhstrt) if none is specified
      end if

   end subroutine ddeabm_initialize