r"""Relaxation experiment.

    The relaxation experiment describes the time-evolution
    of an expectation value starting in a non-equilibrium
    situation.

    Parameters
    ----------
    P : (M, M) ndarray
        Transition matrix
    p0 : (M,) ndarray (optional)
        Initial distribution for a relaxation experiment
    obs : (M,) ndarray
        Observable, represented as vector on state space
    times : list of int (optional)
        List of times at which to compute expectation

    Returns
    -------
    res : ndarray
        Array of expectation value at given times

def relaxation_matvec(P, p0, obs, times=[1]):
    r"""Relaxation experiment.

    The relaxation experiment describes the time-evolution
    of an expectation value starting in a non-equilibrium
    situation.

    Parameters
    ----------
    P : (M, M) ndarray
        Transition matrix
    p0 : (M,) ndarray (optional)
        Initial distribution for a relaxation experiment
    obs : (M,) ndarray
        Observable, represented as vector on state space
    times : list of int (optional)
        List of times at which to compute expectation

    Returns
    -------
    res : ndarray
        Array of expectation value at given times

    """
    times = np.asarray(times)
    """Sort in increasing order"""
    ind = np.argsort(times)
    times = times[ind]

    if times[0] < 0:
        raise ValueError("Times can not be negative")
    dt = times[1:] - times[0:-1]

    nt = len(times)

    relaxations = np.zeros(nt)

    """Propagate obs to initial time"""
    obs_t = 1.0 * obs
    obs_t = propagate(P, obs_t, times[0])
    relaxations[0] = np.dot(p0, obs_t)
    for i in range(nt - 1):
        obs_t = propagate(P, obs_t, dt[i])
        relaxations[i + 1] = np.dot(p0, obs_t)

    """Cast back to original order of time points"""
    relaxations = relaxations[ind]

    return relaxations