The convergence of the discretized version $\max_{x \in \{0, \ldots, N\}} |W(x/N)|$ of $M:=\max_{0 \le t \le 1} |W(t)|$ to $M$ will be very slow -- at the rate of $1/\sqrt N$, according to Korolyuk 1961 and Nagaev 1970 (Korolyuk 1961 apparently exists only in Russian, but can be rather easily read using e.g. Google Translate), 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, using the inverse transform sampling method, 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$:
