def _backward_diff(j, yj, hj, lam):
    """
    Computes through the backward differentiation formula different order
    derivatives of y for a given extrapolation step (for every T{1,j}). In
    other words, it calculates VI.5.43 ref II (step 1) for all kappa for a
    fixed j (for a fixed extrapolation step).

    Used by _interpolate_nonsym.

    Parameters
    ----------
    j : int
        which extrapolation inner-stage is used
    yj : 2D array
        array with 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).
    lam : int
        either 0 or 1, check definition and use in ref II pg 439

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

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

    return rj