Unfortunately the question is in no sense precise. Here is one solution which may be what you want. Assume that a random vector $X = (X_1,\ldots,X_n)$ has the multivariate normal distribution $\mathcal{N}(\mu,\Sigma)$ with mean $\mu = 0$ (the $n$-sphere is centered) and covariance matrix $\Sigma$. The problem seems to be the calculation of
$$\mathbb{P}(X_1^2 + \ldots + X_n^2 \leq \theta_1^2 + \ldots + \theta_n^2)?$$
To simplify things we assume that $\Sigma$ is positive definite (not only semidefinite). Then $\Sigma = A \cdot \Delta \cdot A^T$ with some orthonormal matrix $A$ (i.e. $A \cdot A^T = A^T \cdot A = E_n$, $E_n$ the $n$-dimensional unitmatrix) and a diagonal matrix $\Delta$ with positive diagonal elements $\lambda_1,\ldots,\lambda_n$. For the following you only need $\lambda_i$, $i = 1,\ldots,n$. In R you get $\lambda_i$ with 

    lambda = eigen(Sigma)$values

By assumption $\Sigma = \mathbb{E}XX^T$, hence $\mathbb{E}A^TX(A^TX)^T = A^T\mathbb{E}(XX^T)A = A^T \Sigma A = \Delta$. It follows that $Y = A^TX$ has the covariance matrix $\Delta$. Since the $n$-sphere is not changed by the transformation $A^T$ (here it is important that it is a ball and not an ellipsoid) we can replace the original vector $X$ by the vector $Y = (Y_1,\ldots,Y_n)$ with independent $\mathcal{N}0,\lambda_i)$-distributed $Y_i$ and the problem now is: What is the probability
$$\mathbb{P}(X_1^2 + \ldots + X_n^2 \leq \theta_1^2 + \ldots + \theta_n^2) =\mathbb{P}(Y_1^2 + \ldots + Y_n^2 \leq \theta_1^2 + \ldots + \theta_n^2)?$$
Now $Y_i^2/\lambda_i$ has the distribution $\chi_1^2 = \Gamma_{1/2,1/2}$, (the first parameter is the shape p., the second the rate parameter), hence $Y_i^2$ the distribution $\Gamma_{1/2,\lambda_i/2}$. The distribution function of $Y_1^2 + \ldots + Y_n^2$  can be calculated with the R-library `coga` as
$$\mathbb{P}(Y_1^2 + \ldots + Y_n^2 \leq t) = \text{pcoga}(t,c(1/2,\ldots,1/2),c(\lambda_1/2,\ldots,\lambda_n/2))$$
Inserting $t = \theta_1^2 + \ldots + \theta_n^2$ you get the value you want. There is also an approximate version of this routine: `pcoga_approx`.

To apply coga, you have to open R and then run

`install.library('coda')`

`library(coda)`

`Sigma = matrix(c(1,0.5,0.3,0.5,1,1.4,0.3,1.4,5),nr=3) # only an example`

`lambda = eigen(Sigma)$values`

`n = 3 # dimension of Sigma`

`t = 1.4 # the value you are interested in`

`P = pcoga(t,rep(1/2,n),lambda)`

Then R returns $p = 0.5253459$.