If you just care about the asymptotics, it is indeed just $(1+o(1))\sqrt{\pi n/2}$, where the $o(1)$ term decays like $\tilde{O}(1/n^{1/4})$; this can be done using the approach that I suggested of using the natural embedding to compare to Brownian motion at the obvious stopping times, the fact that the Brownian motion does not fluctuate by more than 1 between stopping times, and the fact that $\vert T_n - n\vert=\tilde{O}(\sqrt{n})$ with very high probability by Bernstein's inequality applied to the (subexponential) increments in stopping times. I can make this more rigorous if needed.
If you really do care about the best constant $C^*$ such that $\mathbb{E}[\max_{k\leq n} S_k]\leq C\sqrt{n}$ for every $n$, not just for large enough $n$, there's a couple of approaches; it's true that $C^*$ is at most twice the optimal constant when you take the absolute values outside the max because $\max_{k\leq n}\vert \sum_{i=1}^k \varepsilon_i\vert\leq \vert\max_{k\leq n} \sum_{i=1}^k \varepsilon_i\vert+\vert\min_{k\leq n} \sum_{i=1}^k \varepsilon_i\vert$, which can be treated exactly via the reflection principle. This should (at least morally) get you something like $2\sqrt{2/\pi}\approx 1.596$. Just for variety, here's a slightly different approach that gives a slightly better constant of $\pi/2\approx 1.571$.
Simply note that if $g_1,\ldots,g_n$ are i.i.d. standard Gaussians, then by Jensen's inequality,
\begin{align*}
\mathbb{E}\left[\max_{k\leq n} \left\vert \sum_{i=1}^k \varepsilon_i \right\vert\right] &= \sqrt{\frac{\pi}{2}} \mathbb{E}\left[\max_{k\leq n} \left\vert \sum_{i=1}^k \varepsilon_i \mathbb{E}[\vert g_i\vert] \right\vert\right]\\
&\leq \sqrt{\frac{\pi}{2}}\mathbb{E}\left[\max_{k\leq n} \left\vert \sum_{i=1}^k \varepsilon_i\vert g_i\vert \right\vert\right]\\
&=\sqrt{\frac{\pi}{2}}\mathbb{E}\left[\max_{k\leq n} \left\vert \sum_{i=1}^k g_i \right\vert\right]\\
&\leq \sqrt{\frac{\pi}{2}}\mathbb{E}\left[\max_{t\leq n} \left\vert B_t\right\vert\right]\\
&=\pi/2\cdot \sqrt{n},
\end{align*}
using the symmetry of standard Gaussians and upper bounding by a coupled Brownian motion and appealing to the MSE solution.