def _compute_H(self, t, index, t2, index2, update_derivatives=False, stationary=False):
if stationary:<EOL><INDENT>raise NotImplementedError("<STR_LIT>")<EOL><DEDENT>Decay = self.decay[index]<EOL>Decay2 = self.decay[index2]<EOL>t_mat = t[:, None]<EOL>t2_mat = t2[None, :]<EOL>if self.delay is not None:<EOL><INDENT>Delay = self.delay[index]<EOL>Delay2 = self.delay[index2]<EOL>t_mat-=Delay[:, None]<EOL>t2_mat-=Delay2[None, :]<EOL><DEDENT>diff_t = (t_mat - t2_mat)<EOL>inv_sigma_diff_t = <NUM_LIT:1.>/self.sigma*diff_t<EOL>half_sigma_decay_i = <NUM_LIT:0.5>*self.sigma*Decay[:, None]<EOL>ln_part_1, sign1 = ln_diff_erfs(half_sigma_decay_i + t2_mat/self.sigma, <EOL>half_sigma_decay_i - inv_sigma_diff_t,<EOL>return_sign=True)<EOL>ln_part_2, sign2 = ln_diff_erfs(half_sigma_decay_i,<EOL>half_sigma_decay_i - t_mat/self.sigma,<EOL>return_sign=True)<EOL>h = sign1*np.exp(half_sigma_decay_i<EOL>*half_sigma_decay_i<EOL>-Decay[:, None]*diff_t+ln_part_1<EOL>-np.log(Decay[:, None] + Decay2[None, :]))<EOL>h -= sign2*np.exp(half_sigma_decay_i*half_sigma_decay_i<EOL>-Decay[:, None]*t_mat-Decay2[None, :]*t2_mat+ln_part_2<EOL>-np.log(Decay[:, None] + Decay2[None, :]))<EOL>if update_derivatives:<EOL><INDENT>sigma2 = self.sigma*self.sigma<EOL>dh_ddecay = ((<NUM_LIT:0.5>*Decay[:, None]*sigma2*(Decay[:, None] + Decay2[None, :])-<NUM_LIT:1>)*h<EOL>+ (-diff_t*sign1*np.exp(<EOL>half_sigma_decay_i*half_sigma_decay_i-Decay[:, None]*diff_t+ln_part_1<EOL>)<EOL>+t_mat*sign2*np.exp(<EOL>half_sigma_decay_i*half_sigma_decay_i-Decay[:, None]*t_mat<EOL>- Decay2*t2_mat+ln_part_2))<EOL>+self.sigma/np.sqrt(np.pi)*(<EOL>-np.exp(<EOL>-diff_t*diff_t/sigma2<EOL>)+np.exp(<EOL>-t2_mat*t2_mat/sigma2-Decay[:, None]*t_mat<EOL>)+np.exp(<EOL>-t_mat*t_mat/sigma2-Decay2[None, :]*t2_mat<EOL>)-np.exp(<EOL>-(Decay[:, None]*t_mat + Decay2[None, :]*t2_mat)<EOL>)<EOL>))<EOL>self._dh_ddecay = (dh_ddecay/(Decay[:, None]+Decay2[None, :])).real<EOL>dh_ddecay2 = (t2_mat*sign2<EOL>*np.exp(<EOL>half_sigma_decay_i*half_sigma_decay_i<EOL>-(Decay[:, None]*t_mat + Decay2[None, :]*t2_mat)<EOL>+ln_part_2<EOL>)<EOL>-h)<EOL>self._dh_ddecay2 = (dh_ddecay/(Decay[:, None] + Decay2[None, :])).real<EOL>self._dh_dsigma = (half_sigma_decay_i*Decay[:, None]*h<EOL>+ <NUM_LIT:2>/(np.sqrt(np.pi)<EOL>*(Decay[:, None]+Decay2[None, :]))<EOL>*((-diff_t/sigma2-Decay[:, None]/<NUM_LIT:2>)<EOL>*np.exp(-diff_t*diff_t/sigma2)<EOL>+ (-t2_mat/sigma2+Decay[:, None]/<NUM_LIT:2>)<EOL>*np.exp(-t2_mat*t2_mat/sigma2-Decay[:, None]*t_mat) <EOL>- (-t_mat/sigma2-Decay[:, None]/<NUM_LIT:2>) <EOL>*np.exp(-t_mat*t_mat/sigma2-Decay2[None, :]*t2_mat) <EOL>- Decay[:, None]/<NUM_LIT:2><EOL>*np.exp(-(Decay[:, None]*t_mat+Decay2[None, :]*t2_mat))))<EOL><DEDENT>return h<EOL>
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