def _compute_H(self, t, index, t2, index2, update_derivatives=False, stationary=False):
        """Helper function for computing part of the ode1 covariance function.

        :param t: first time input.
        :type t: array
        :param index: Indices of first output.
        :type index: array of int
        :param t2: second time input.
        :type t2: array
        :param index2: Indices of second output.
        :type index2: array of int
        :param update_derivatives: whether to update derivatives (default is False)
        :return h : result of this subcomponent of the kernel for the given values.
        :rtype: ndarray
"""

        if stationary:
            raise NotImplementedError, "Error, stationary version of this covariance not yet implemented."
        # Vector of decays and delays associated with each output.
        Decay = self.decay[index]
        Decay2 = self.decay[index2]
        t_mat = t[:, None]
        t2_mat = t2[None, :]
        if self.delay is not None:
            Delay = self.delay[index]
            Delay2 = self.delay[index2]
            t_mat-=Delay[:, None]
            t2_mat-=Delay2[None, :]

        diff_t = (t_mat - t2_mat)
        inv_sigma_diff_t = 1./self.sigma*diff_t
        half_sigma_decay_i = 0.5*self.sigma*Decay[:, None]

        ln_part_1, sign1 = ln_diff_erfs(half_sigma_decay_i + t2_mat/self.sigma, 
                                        half_sigma_decay_i - inv_sigma_diff_t,
                                        return_sign=True)
        ln_part_2, sign2 = ln_diff_erfs(half_sigma_decay_i,
                                        half_sigma_decay_i - t_mat/self.sigma,
                                        return_sign=True)

        h = sign1*np.exp(half_sigma_decay_i
                         *half_sigma_decay_i
                         -Decay[:, None]*diff_t+ln_part_1
                         -np.log(Decay[:, None] + Decay2[None, :]))
        h -= sign2*np.exp(half_sigma_decay_i*half_sigma_decay_i
                          -Decay[:, None]*t_mat-Decay2[None, :]*t2_mat+ln_part_2
                          -np.log(Decay[:, None] + Decay2[None, :]))

        if update_derivatives:
            sigma2 = self.sigma*self.sigma
            # Update ith decay gradient

            dh_ddecay = ((0.5*Decay[:, None]*sigma2*(Decay[:, None] + Decay2[None, :])-1)*h
                         + (-diff_t*sign1*np.exp(
                half_sigma_decay_i*half_sigma_decay_i-Decay[:, None]*diff_t+ln_part_1
                )
                            +t_mat*sign2*np.exp(
                half_sigma_decay_i*half_sigma_decay_i-Decay[:, None]*t_mat
                - Decay2*t2_mat+ln_part_2))
                         +self.sigma/np.sqrt(np.pi)*(
                -np.exp(
                -diff_t*diff_t/sigma2
                )+np.exp(
                -t2_mat*t2_mat/sigma2-Decay[:, None]*t_mat
                )+np.exp(
                -t_mat*t_mat/sigma2-Decay2[None, :]*t2_mat
                )-np.exp(
                -(Decay[:, None]*t_mat + Decay2[None, :]*t2_mat)
                )
                ))
            self._dh_ddecay = (dh_ddecay/(Decay[:, None]+Decay2[None, :])).real
            
            # Update jth decay gradient
            dh_ddecay2 = (t2_mat*sign2
                         *np.exp(
                half_sigma_decay_i*half_sigma_decay_i
                -(Decay[:, None]*t_mat + Decay2[None, :]*t2_mat)
                +ln_part_2
                )
                         -h)
            self._dh_ddecay2 = (dh_ddecay/(Decay[:, None] + Decay2[None, :])).real
            
            # Update sigma gradient
            self._dh_dsigma = (half_sigma_decay_i*Decay[:, None]*h
                               + 2/(np.sqrt(np.pi)
                                    *(Decay[:, None]+Decay2[None, :]))
                               *((-diff_t/sigma2-Decay[:, None]/2)
                                 *np.exp(-diff_t*diff_t/sigma2)
                                 + (-t2_mat/sigma2+Decay[:, None]/2)
                                 *np.exp(-t2_mat*t2_mat/sigma2-Decay[:, None]*t_mat) 
                                 - (-t_mat/sigma2-Decay[:, None]/2) 
                                 *np.exp(-t_mat*t_mat/sigma2-Decay2[None, :]*t2_mat) 
                                 - Decay[:, None]/2
                                 *np.exp(-(Decay[:, None]*t_mat+Decay2[None, :]*t2_mat))))
                
        return h