<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def correlation_matvec(P, obs1, obs2=None, times=[1]):
 r"""Time-correlation for equilibrium experiment - via matrix vector products. Parameters P : (M, M) ndarray Transition matrix obs1 : (M,) ndarray Observable, represented as vector on state space obs2 : (M,) ndarray (optional) Second observable, for cross-correlations times : list of int (optional) List of times (in tau) at which to compute correlation Returns ------- correlations : ndarray Correlation values at given times """

    if obs2 is None:
        obs2 = obs1

    """Compute stationary vector"""
    mu = statdist(P)
    obs1mu = mu * obs1

    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)

    correlations = np.zeros(nt)

    """Propagate obs2 to initial time"""
    obs2_t = 1.0 * obs2
    obs2_t = propagate(P, obs2_t, times[0])
    correlations[0] = np.dot(obs1mu, obs2_t)
    for i in range(nt - 1):
        obs2_t = propagate(P, obs2_t, dt[i])
        correlations[i + 1] = np.dot(obs1mu, obs2_t)

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

    return correlations