Probability that one RV will exceed many others Assume the $1 \times N$ vector
$\mathbf X = [X_1, X_2, \ldots , X_N]$
contains i.i.d. normal samples such that $\mathbf X$ has a multivariate normal distribution. Now assume another random variable $Y$ also has a normal distribution, independent from the samples in $\mathbf X$, such that $Y$ has a mean and variance not necessarily equal to that of $\mathbf X$.
What is the probability that $Y$ will exceed all the values in $\mathbf X$ concurrently (or jointly)? In other words what is the probability that the value of $Y$ will be the maximum value of $[X_1, X_2, \ldots , X_N, Y]$
Possible approach?
In the simplest case $N=1$, the problem reduces to 
$p(Y>X_1)$ = $p(Y-X_1>0)$
which can be calculated by convolution of $p(Y)$ and $p(-X_1)$ and integrating the resulting PDF from $0$ to $\infty$. When $Y$ is also i.i.d. w.r.t the samples of $\mathbf X$, the resulting PDF is zero-mean, such that $p(Y>X_1)$ = 0.5. This makes intuitively sense, since the probability that one RV will exceed another i.i.d. RV is 0.5. The probability
$p(Y>X_1, Y>X_2, \ldots, Y>X_N )$
should therefore equal $1/(N+1)$ if $Y$ is i.i.d. w.r.t. $\mathbf X$. When $N>1$, the problem does not seem straightforward (maybe it is?).
 A: Assume without loss of generality that each $X_i$ is standard normal and that $Y$ is normal with mean $\mu$ and variance $\sigma^2$. By definition, for every $y$,
$$
P[y\gt\max(X_1,\ldots,X_N)]=\Phi(y)^N,
$$
hence, by the independence of $Y$ from $(X_i)_i$,
$$
P[Y\gt\max(X_1,\ldots,X_N)]=E[\Phi(Y)^N]=\int_\mathbb R\Phi(\mu+\sigma x)^N\varphi(x)\mathrm dx,
$$
where
$$
\varphi(x)=\frac1{\sqrt{2\pi}}\mathrm e^{-x^2/2},\qquad\Phi(x)=\int_{-\infty}^x\varphi(z)\mathrm dz.
$$
Except when $\mu=0$, $\sigma^2=1$, I see no simplification. An equivalent formula is
$$
P[Y\gt\max(X_1,\ldots,X_N)]=\int_\mathbb R\Phi(x)^N\varphi\left(\frac{x-\mu}\sigma\right)\frac{\mathrm dx}\sigma.
$$
A: The case $n=1$ can be done more explicitly than you say. If $X \sim \mathcal N(\mu_X,\sigma_X^2), Y \sim \mathcal N(\mu_Y,\sigma_Y^2)$ are independent then $Y-X \sim \mathcal N(\mu_Y-\mu_X, \sigma_X^2 + \sigma_Y^2)$ so $P(Y - X \gt 0) = 1 -\Phi(\frac{\mu_X - \mu_Y}{\sigma_X^2 + \sigma_Y^2}) = \Phi(\frac{\mu_Y-\mu_X}{\sigma_X^2+\sigma_Y^2}).$
If you are interested in a particular value of $n$, then I think you should use numerical integration which works quite well. You may also be able to get some asymptotics. A lot is known about the distribution of order statistics including the maximum of $n$ IID normal distributions. For example, the mean value for the maximum of $n$ standard normals is approximately $\Phi^{-1}(\frac{n-0.375}{n+0.25}).$ See also the accepted answer to "Expectation of the maximum of gaussian random variables" and the reference in the comments.
The maximum of $n$ IID normals is not normal. However, it may be ok to approximate it as a normal distribution, and use the case $n=1$, since it doesn't start too far away from a normal distribution (it is continuous and unimodal) and you are essentially convolving it with one Gaussian afterwards. This approximation should work better closer to the mean.
A: The probability that $Y$ will exceed all the values in $\mathbf X$ can be written as (see answer of Didier Piau above)
$p_m = p(Y > {\rm max}(\mathbf X)) = p(Y - {\rm max}(\mathbf X) > 0)$
The density of ${\rm max}(\mathbf X)$ can be obtained using (Papoulis p. 193)
$f_{\rm max}(z) = N f_X(z) F_X^{N-1}(z)$
with $f_X$ the density of $X_i$ within $\mathbf X$ and $F_X$ the associated distribution function (CDF).
The density of the difference between the RVs $Y$ and ${\rm max}(\mathbf X)$ can be obtained using (Papoulis p. 185)
$p(Y - {\rm max}(\mathbf X)) = f_z(z) = \int_{-\infty}^{\infty} f_Y(z+w) f_{max}(w)~dw$
which is simply the convolution between $f_Y = p(Y)$ and $f_{max} = p({\rm max}(\mathbf X))$. The final answer is then the area under $f_z(z)$ beyond $0$ which can be written
$p_m = 1-F_z(0)$
