Numerically computing Gaussian integral

Question

Given a trivariate Gaussian random variable

$$X\sim N(\mu,\Sigma)$$

I'm trying to compute the expected values

$$\text{E}\Big[\frac{e^{x_i}}{e^{x_1}+e^{x_2}+e^{x_3}}\Big]=\frac{1}{\sqrt{(2\pi)^3\det\Sigma}}\int_{{\bf R}^3}\frac{e^{x_i}}{e^{x_1}+e^{x_2}+e^{x_3}}e^{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)}$$

for $i=1,2,3$.

I've used several different cubature packages to numerically compute this integral, but for certain values of $\Sigma$ it either fails to converge or it returns garbage output.

Here is an example of a parameter set which causes problems:

\begin{align} \mu&=(-3.4757907, 4.7485485, 0.7119134)\notag\\\\ \Sigma&=\begin{bmatrix} 0.0336 & -0.0215 & -0.0130 \\ -0.0215 & 0.209 & -0.0922 \\ -0.0130 & -0.0922 & 0.1558 \end{bmatrix}\notag \end{align}

It's not entirely clear why the numerical integration packages are failing, I assumed it was either underflow of the denominator or overflow in the exponentials, but when I set $\Sigma=I$ I'm able to successfully perform the integration, which makes me think it might be something else.

I considered forming the approximation

$$e^{x_1}+e^{x_2}+e^{x_3}\approx e^{\max(x)},$$

so that I could push all the terms into a single exponential, but it didn't help.

Obviously if it comes to it I can compute these expected values using Monte Carlo methods, however I'd prefer to avoid that if possible.

using Cubature

B = 200

mu = [-3.4757907, 4.7485485, 0.7119134]
sigma = [0.0336 -0.0215 -0.0130; -0.0215 0.209 -0.0922; -0.0130 -0.0922 0.1558]

inv_sigma = inv(sigma)
C = sqrt(det(2pi * sigma))

function func1(x)
    a = -0.5 * (x - mu)' * inv_sigma * (x - mu)
    f = exp(x[1] + a) / (C * sum(exp.(x)))
    return(f)
end

function func2(x)
    a = -0.5 * (x - mu)' * inv_sigma * (x - mu)
    f = exp(x[2] + a) / (C * sum(exp.(x)))
    return(f)
end

function func3(x)
    a = -0.5 * (x - mu)' * inv_sigma * (x - mu)
    f = exp(x[3] + a) / (C * sum(exp.(x)))
    return(f)
end

p1 = hcubature(func1, [-B,-B,-B], [B,B,B])[1]
p2 = hcubature(func2, [-B,-B,-B], [B,B,B])[1]
p3 = hcubature(func3, [-B,-B,-B], [B,B,B])[1]