I think you can start with the process of uniformly selecting a point on (n − 1)-sphere and then consider the original samples as independent.

So we draw a vector $Y$ with $n$ coordinates taken i.i.d. from the normal distribution $Y_i \sim \mathcal{N}(0, 1)$. The CDF of the absolute maximum of those is given by taking the maximum of $n$ points from the half normal distribution  $F_{||Y||_{\infty}}(x)={\operatorname{erf}\left( \frac{ x }{\sqrt 2} \right) }^{n}$. We do need to normalize the points, so they'll be on the unit sphere. Given that the positive squared root of sum of squared of normal variables has a Chi distribution with $n$ degress of freedom, $||Y||_2\sim \chi_n$, we have:

$$P(||X||_{\infty} < r) = P(||Y||_{\infty}/|Y||_2 < r) = \\\int_{0}^{\infty} P(||Y||_{\infty}/ ||Y||_2 < r\; |\; ||Y||_2 = x)  f_{\chi_n}(x)\;dx = \\ 
\\\int_{0}^{\infty} F_{||Y||_{\infty}}(x \cdot r) f_{\chi_n}(x)\;dx = \\\int_{0}^{\infty} {\operatorname{erf}\left( \frac{ x \cdot r}{\sqrt 2} \right) }^{n}  \cdot \dfrac{x^{n-1}e^{-x^2/2}}{2^{n/2-1}\Gamma\left(\frac{n}{2}\right)}\; dx$$