Running a simulation, as Tom K suggests, is one option, but the computation becomes prohibitively time consuming if you have to repeat the simulation with different parameters. Here is a faster option:
Let F1 denote the marginal cdf of X1, and similarly for X2,...,Xn. Let f1 denote the marginal pdf of X1, and similarly for X2,...,Xn. Finally, let f denote the joint pdf of X1,...,Xn, which by independence is equal to the product of the marginal pdfs.
The desired probability is
P(X1 > max{X2,...,Xn})
= P(X1 > X2,...,X1>Xn)
=∫(∫···∫f(x1,x2,...,xn)dx2···dxn)dx1
=∫(∫···∫f(x1)f(x2)···f(xn)dx2···dxn)dx1
where the outer-most integral is from -∞ and +∞, while all other integrals are from -∞ to x1. Rearranging factors and noting that Fi(x1) = ∫-∞x1 fi(xi)dxi for each i=2,...,n, we get
P(X1 > max{X2,...,Xn})
= ∫f1(x1)F2(x1)···Fn(x1)dx1
with limits -∞ and +∞. Notice that the first factor is a pdf, while all of the other factors are cdfs.
This result applies to any sequence of independent random variables for which the pdfs/cdfs are known. A good numerical integrator, such as in Mathematica or Maple, can approximate the probability very quickly. I wouldn't recommend R, however. Its integrate function can suffer from very poor precision, including in the special case you mentioned. Wolframalpha works just fine if you don't have Mathematica.
For demonstration, let's assume X1~N(40,10), X2~N(36,8), X3~N(45,12), X4~N(35,15).
Go to the wolframalpha website and query
PDF[NormalDistribution[40,10], x] * CDF[NormalDistribution[36,8], x] * CDF[NormalDistribution[45,12], x] * CDF[NormalDistribution[35,15], x]
This will give the function to be integrated. Scroll down, hover the cursor over "Exact result", click on "Plain text" and then copy the Wolfram Language plain text output. Next, query
Integrate( ____ ,x,-inf,inf)
but in the blank space, paste the text you copied. Wolfram should respond in a few seconds with the approximate probability 0.24076.
