$\newcommand\la\lambda$By rescaling, it is enough to consider the case $\la=1$ (multiplying the resulting expectation by $\la$ at the end). By symmetry, $Y:=\max X$ and $-\min X$ are identically distributed. So, for $M_n:=M$ we have
\begin{equation*}
EM_n=2EY. \tag{10}\label{10}
\end{equation*}
Next, $Y=\max(0,Y)-\max(0,-Y)$ and hence
\begin{equation*}
EY=I_1-I_2, \tag{20}\label{20}
\end{equation*}
where
\begin{equation*}
\begin{aligned}
I_1&:=E\max(0,Y)=\int_0^\infty dx\,P(Y>x)=\int_0^\infty dx\,(1-(1-F(x)^n),\\
I_2&:=E\max(0,-Y)=\int_0^\infty dx\,P(-Y>x)=\int_0^\infty dx\,F(-x)^n.
\end{aligned}
\end{equation*}
Here I use a standard trick:
\begin{equation*}
\int_0^\infty dx\,P(Y>x)=\int_0^\infty dx\,E\,1(Y>x) \\
=E\int_0^\infty dx\,1(Y>x)
=E\int_0^{\max(0,Y)} dx=E\max(0,Y).
\end{equation*}
Further,
\begin{equation*}
I_2=\int_0^\infty dx\,(e^{-x}/2)^n=\frac{2^{-n}}n \tag{30}\label{30}
\end{equation*}
and
\begin{equation*}
\begin{aligned}
I_1&=\int_0^\infty dx\,(1-(1-e^{-x}/2)^n) \\
&=-\int_0^\infty dx\,\sum_{j=1}^n\binom nj(-e^{-x}/2)^j \\
&=-\sum_{j=1}^n\binom nj\frac{(-1/2)^j}j
=\frac{n}{2} \, _3F_2\left(1,1,1-n;2,2;\frac{1}{2}\right),
\end{aligned}
\end{equation*}
where $_3F_2$ is the hypergeometric function.
Thus, for $\la=1$,
\begin{equation*}
EM_n=n \, _3F_2\Big(1,1,1-n;2,2;\frac{1}{2}\Big)-\frac1{n2^{n-1}}
\end{equation*}
and
\begin{equation*}
EM_n=\la\Big(n \, _3F_2\Big(1,1,1-n;2,2;\frac{1}{2}\Big)-\frac1{n2^{n-1}}\Big)
\end{equation*}
for any real $\la>0$.
It follows from the previous answer (thanks to Matt F. for reminding us of it), formulas \eqref{10}, \eqref{20}, \eqref{30}, and the "rescaling" remark that
\begin{equation*}
l_n:=\gamma+\ln(n/2)-\frac1{n{2^n}}<\frac{EM_n}{2\la}<u_n:=1+\ln(n/2)-\frac1{n{2^n}},
\end{equation*}
where $\gamma=0.577\ldots$ is the Euler gamma.
It is clear now that
\begin{equation*}
\frac{EM_n}{2\la}\sim l_n\sim u_n\sim\ln(n/2)\sim\ln n,
\end{equation*}
so that the upper bound $u_n$ and the lower bound $l_n$ on $\frac{EM_n}{2\la}$ are each asymptotically exact.
