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)

    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.

    real(wp),parameter :: default_tol = 100*epsilon(1.0_wp) !! if tols not specified

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

    ! 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

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

    end subroutine initialize_variable_step