<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _matrix_sigma_eps(self, donor_catchments):
 """ Return sampling error coveriance matrix Sigma eta Methodology source: Kjeldsen & Jones 2009, eq 9 :param donor_catchments: Catchments to use as donors :type donor_catchments: list of :class:`Catchment` :return: 2-Dimensional, symmetric covariance matrix :rtype: :class:`numpy.ndarray` """

        p = len(donor_catchments)
        sigma = np.empty((p, p))
        for i in range(p):
            beta_i = self._beta(donor_catchments[i])
            n_i = donor_catchments[i].amax_records_end() - donor_catchments[i].amax_records_start() + 1
            for j in range(p):
                beta_j = self._beta(donor_catchments[j])
                n_j = donor_catchments[j].amax_records_end() - donor_catchments[j].amax_records_start() + 1
                rho_ij = self._lnqmed_corr(donor_catchments[i], donor_catchments[j])
                n_ij = min(donor_catchments[i].amax_records_end(), donor_catchments[j].amax_records_end()) - \
                       max(donor_catchments[i].amax_records_start(), donor_catchments[j].amax_records_start()) + 1
                sigma[i, j] = 4 * beta_i * beta_j * n_ij / n_i / n_j * rho_ij
        return sigma