The way you propose to simulate $M:=\max_{0 \le t \le 1} |W(t)|$ is very inefficient. The convergence of the discretized version (say $M_N$) of $M$ to $M$ will be very slow, probably at the rate of $1/\sqrt N$ or so, and then you will have to simulate a very large number $N$ of normal random variables -- just to get one realization of $M$.
A much more efficient way is to use the following explicit expression for the cdf of $M$: \begin{equation} P(M\le x)=\frac4\pi\, \sum _{k=0}^{\infty } \frac{(-1)^k}{2 k+1}\, \exp \left(-\frac{(2 k+1)^2 \pi ^2}{8 x^2}\right) \tag{1} \end{equation} for $x>0$; see e.g. page 3. The series in (1) converges very fast unless $x$ is very large. On the other hand, if $x$ is large, then $P(M\le x)$ is very close to $1$.
More specifically, if $x>x_0$, where $x_0=5.8471\ldots$ is the root of the equation $4(1-\Phi(x_0))=1/10^8$ and $\Phi$ is the standard normal cdf, then
\begin{equation}
0<1-P(M\le x)=P(M>x)\le P(M^+>x)+P(M^->x)=2P(M^+>x)=4(1-\Phi(x))<1/10^8,
\end{equation}
where $M^\pm:=\max_{0\le t\le1}(\pm W(t))$; the last displayed equality follows immediately
from another formula on the same page 3.
On the other hand, for $x\in(0,x_0]$, retaining only the first $11$ summands of the alternating series in (1), we will get the approximation \begin{equation} P(M\le x)\approx F(x):=\frac4\pi\, \sum _{k=0}^{10} \frac{(-1)^k}{2 k+1}\, \exp \left(-\frac{(2 k+1)^2 \pi ^2}{8 x^2}\right) \end{equation} with an error less in absolute value than \begin{equation} \frac4\pi\,\frac1{2\times11+1}\, \exp \Big(-\frac{(2\times11+1)^2 \pi ^2}{8x_0^2}\Big)<3/10^{10}. \end{equation}
Now we can get very close approximations $m_1,m_2,\dots$ to iid realizations of $M$ as solutions of the equations $F(m_i)=u_i$ for $m_i$, where $u_1,u_2,\dots$ are iid realizations of a random variable uniformy distributed between $0$ and $1$. It then takes under a second in Mathematica to simulate $1000$ very close approximations $m_1,\dots,m_{1000}$ to iid realizations of $M$:
