def _centered_diff(j, f_yj, hj):
    """
    Compute approximate derivatives of y using centered differences.
    Specifically, calculate II.9.39 of Ref. I (step 1) for all kappa for a
    fixed j (i.e. for a fixed extrapolation step).

    Used by _interpolate_sym.

    Parameters
    ----------
    j : int
        which extrapolation inner stage is used
    f_yj : 2D array
        array with all the function evaluations of all the intermediate
        solution values obtained to calculate each T_{j,1}.
    hj : array
        inner step taken in the j-th extrapolation step, H/nj (II.9.1 ref I).

    dj : 2D array
        contains, for each kappa=1...2j, the kappa-th derivative of y estimated
        using the j-th extrapolation step values
    """
    max_order = 2*j
    nj = len(f_yj) - 1
    coeff = [1,1]
    dj = (max_order+1)*[None]
    dj[1] = 1*f_yj[nj//2]  # First derivative
    dj[2] = (f_yj[nj//2+1] - f_yj[nj//2-1])/(2*hj)  # Second derivative

    for order in range(2,max_order):  # Higher-order derivatives
        # Binomial coefficients
        coeff = [1] + [coeff[jj] + coeff[jj+1] for jj in range(len(coeff)-1)] + [1]
        index = range(nj//2+order, nj//2-order-2, -2)
        sum_ = sum([((-1)**i)*coeff[i]*f_yj[index[i]] for i in range(order+1)])
        dj[order+1] = sum_ / (2*hj)**order

    return dj