I found here that
The cdf for the max is the cdf for the normal raised to the power of the sample size. [say, $n$]
Since the cdf for the normal is the error function,
$$F_X(x)=\frac{1}{\sigma\sqrt{ 2\pi}} \int_{-\infty}^x \exp{\left(\frac{-(t -\mu)^2}{2\sigma^2} \right)}\mathrm dt=\frac{1}{2}\Bigg[1+ \mathrm{erf} \left( \frac{x-\mu}{\sigma \sqrt 2}\right)\Bigg],$$
and the pdf of the standard normal is
$$f(x) = \frac{e^{-x^{2}/2}} {\sqrt{2\pi}}$$
I presume there is no other workaround to
$$f_{\max(\mathrm X)}(x_i)=n \quad \frac{e^{-x^{2}/2}} {\sqrt{2\pi}}\quad \frac{1}{2}\Bigg[1+ \mathrm{erf} \left( \frac{x-\mu}{\sigma \sqrt 2}\right)\Bigg]^{n-1}$$
Would this be correct? It doesn't get any simpler?