initialize_variable_step Subroutine

    subroutine initialize_variable_step(me,n,f,rtol,atol,stepsize_method,&
                                        hinit_method,report,g,stop_on_errors,&
                                        max_number_of_steps,report_rate,&
                                        solver,dense_output,dense_output_step)

    implicit none

    class(rk_variable_step_class),intent(inout) :: me
    integer,intent(in)                          :: n               !! number of equations
    procedure(deriv_func)                       :: f               !! derivative function
    real(wp),dimension(:),intent(in),optional   :: rtol            !! relative tolerance (if size=1,
                                                                   !! then same tol used for all
                                                                   !! equations). If not present, a default
                                                                   !! of `100*eps` is used
    real(wp),dimension(:),intent(in),optional   :: atol            !! absolute tolerance (if size=1,
                                                                   !! then same tol used for all
                                                                   !! equations). If not present, a default
                                                                   !! of `100*eps` is used
    type(stepsize_class),intent(in),optional    :: stepsize_method !! method for varying the step size
    integer,intent(in),optional                 :: hinit_method    !! which method (1 or 2) to use for
                                                                   !! automatic initial step size
                                                                   !! computation.
                                                                   !! 1 = use `hstart`, 2 = use `hinit`.
    procedure(report_func),optional             :: report          !! for reporting the steps
    procedure(event_func),optional              :: g               !! for stopping at an event
    logical,intent(in),optional :: stop_on_errors !! stop the program for
                                                  !! any errors (default is False)
    integer,intent(in),optional :: max_number_of_steps !! max number of steps allowed
    integer,intent(in),optional :: report_rate !! how often to call report function:
                                               !! `0` : no reporting (same as not associating `report`),
                                               !! `1` : report every point,
                                               !! `2` : report every other point, etc.
                                               !! The first and last point are always reported.
    class(root_method),intent(in),optional :: solver !! the root-finding method to use for even finding.
                                                     !! if not present, then `brent_solver` is used.
    logical,intent(in),optional :: dense_output !! to enable the dense output (default is False, unless
                                                !! `dense_output_step` is present and not zero).
                                                !! If enabled, [[interpolate]] can be used (e.g., from the
                                                !! `report` function), and the dense output is also used for
                                                !! event location in [[integrate_to_event_variable_step]].
                                                !! The method must have a dense output. When not enabled,
                                                !! there is no overhead for the dense output.
    real(wp),intent(in),optional :: dense_output_step !! if present and not zero, the points are reported
                                                      !! at this fixed time interval (starting from the initial
                                                      !! time), using the dense output of the method,
                                                      !! rather than at the steps taken by the integrator.
                                                      !! The first and last point are also reported (so the
                                                      !! last interval will generally be shorter).
                                                      !! This is only used if `report` is present
                                                      !! and `report_rate>0` (`report_rate` is otherwise ignored).
                                                      !! The method must have a dense output. This automatically
                                                      !! enables the dense output (it is an error if
                                                      !! `dense_output=.false.` is also specified).

    real(wp),parameter :: default_tol = 100*epsilon(1.0_wp) !! if tols not specified
    type(rklib_properties) :: props !! to check if the method has a dense output
    logical :: use_dense_step !! if `dense_output_step` was specified
    logical :: enable_dense !! if the dense output is to be enabled

    ! dense output:
    use_dense_step = .false.
    if (present(dense_output_step)) use_dense_step = abs(dense_output_step) > zero
    if (present(dense_output)) then
        enable_dense = dense_output
    else
        enable_dense = use_dense_step
    end if

    ! base init:
    call me%init(n,f,report,g,stop_on_errors,max_number_of_steps,report_rate,solver,enable_dense)

    ! variable-step specific inputs:
    if (allocated(me%rtol)) deallocate(me%rtol)
    if (allocated(me%atol)) deallocate(me%atol)
    allocate(me%rtol(n))
    allocate(me%atol(n))

    if (present(rtol)) then
        if (size(rtol)==1) then
            me%rtol = abs(rtol(1)) !use this for all equations
        else if (size(rtol)==n) then
            me%rtol = abs(rtol)
        else
            call me%raise_exception(RKLIB_ERROR_INVALID_RTOL_SIZE)
        end if
    else
        me%rtol = default_tol
    end if

    if (present(atol)) then
        if (size(atol)==1) then
            me%atol = abs(atol(1)) !use this for all equations
        else if (size(atol)==n) then
            me%atol = abs(atol)
        else
            call me%raise_exception(RKLIB_ERROR_INVALID_ATOL_SIZE)
        end if
    else
        me%atol = default_tol
    end if

    if (present(hinit_method)) then
        if (any(hinit_method == [1,2])) then
            me%hinit_method = hinit_method
        else
            call me%raise_exception(RKLIB_ERROR_INVALID_HINIT_METHOD)
            return
        end if
    end if
    if (present(stepsize_method)) me%stepsize_method = stepsize_method

    ! check the dense output inputs:
    if (use_dense_step .and. .not. enable_dense) then
        call me%raise_exception(RKLIB_ERROR_DENSE_STEP_DISABLED)
        return
    end if
    if (enable_dense) then
        props = me%properties()
        if (.not. props%dense) then
            call me%raise_exception(RKLIB_ERROR_NO_DENSE_OUTPUT)
            return
        end if
    end if
    if (use_dense_step) me%dense_output_step = abs(dense_output_step)

    ! reset internal variables:
    me%num_rejected_steps = 0
    me%last_accepted_step_size = zero

    end subroutine initialize_variable_step