def diff(xi, yi, order=1) -> np.ndarray:
    """Take the numerical derivative of a 1D array.

    Output is mapped onto the original coordinates  using linear interpolation.
    Expects monotonic xi values.

    Parameters
    ----------
    xi : 1D array-like
        Coordinates.
    yi : 1D array-like
        Values.
    order : positive integer (optional)
        Order of differentiation.

    Returns
    -------
    1D numpy array
        Numerical derivative. Has the same shape as the input arrays.
    """
    yi = np.array(yi).copy()
    flip = False
    if xi[-1] < xi[0]:
        xi = np.flipud(xi.copy())
        yi = np.flipud(yi)
        flip = True
    midpoints = (xi[1:] + xi[:-1]) / 2
    for _ in range(order):
        d = np.diff(yi)
        d /= np.diff(xi)
        yi = np.interp(xi, midpoints, d)
    if flip:
        yi = np.flipud(yi)
    return yi
Take the numerical derivative of a 1D array.

    Output is mapped onto the original coordinates  using linear interpolation.
    Expects monotonic xi values.

    Parameters
    ----------
    xi : 1D array-like
        Coordinates.
    yi : 1D array-like
        Values.
    order : positive integer (optional)
        Order of differentiation.

    Returns
    -------
    1D numpy array
        Numerical derivative. Has the same shape as the input arrays.