Covariance matrix of integral of multivariate normal distribution

If $t = [t_0, t_1, dots, t_{N-1}] in mathbb{R}^N$ with $t_i sim N(mu_i, sigma_i^2)$ and its covariance matrix $C in mathbb{R}^{N times N}$ where $C_{ij} = Cov(t_i, t_j)$ is given

If I define a multi bernullian distribution $y = t > 0$ where $y_i = t_i > 0$ (1 if true, zero if false)

$E[y] = 1 – frac{1}{2}(1 +erf(-frac{E[t]}{sqrt{2}Var[t]}))$

Where $Var[t] = [C_{00}, C_{11}, dots, C_{(N-1)(N-1)}]$

How can I compute the covariance matrix of $y$? Did I compute $E[y]$ correctly?