docs: 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