I'm pretty sure this is unknown, though it would be great if I'm wrong.

Asymptotically there is a function $A(n)$ such that $M(n,n/2+t)\sim  e^{-2t^2}A(n)$ for $t=o(n^{1/2})$, which follows from [this paper][1] of Canfield and McKay.

It is also known to be true for $n\le 20$.

  [1]: http://www.combinatorics.org/ojs/index.php/eljc/article/view/v12i1r29

To spell out the asymptotics a bit more, the paper cited above shows that
$$ M(n,k) = (e^{-1/2}+o(1)) \binom {n}{k}^{\!2n} \bigl( \lambda^\lambda(1-\lambda)^{1-\lambda}\bigr)^{n^2}, \qquad(*)$$
as $n\to\infty$, where $\lambda=k/n$ and $cn/\log n\le k\le n-cn/\log n$ for a particular constant $c$. The value $c=1/3$ will do. This expression without the error term is unimodal for $0\le k\le n$ and even with the error term it is unimodal if $n$ is large enough. The same formula holds for $k=o(n^{1/2})$ (McKay and Wang, 2003). A new method of Liebenau and Wormald will (I'm 100% confident) show that $(*)$ is also true for the intermediate ranges of $k$, but it is not published yet. Then we will know that $M(n,k)$ is unimodal in $k$ provided $n$ is large enough.

If we just want to know whether $M(n,k)\le M(n,\lfloor n/2\rfloor)$, as asked, and don't care about unimodality, then what remains asymptotically is to show that
$M(n,k)\lt M(n,\lfloor n/2\rfloor)$ for $k\le \frac13 n/\log n$. This follows from the pairing (configuration) model for random $k$-regular bipartite graphs; namely
$$ M(n,k) \le \frac{ (nk)! }{ (k!)^{2n} }.\qquad(\#)$$
For large $k$, $(\#)$ is a pretty terrible bound, but if I didn't miscalculate it is sufficient to prove that $M(n,k)\lt M(n,\lfloor n/2\rfloor)$ for $k\le \frac13 n/\log n$.

That completes the proof that $M(n,\lfloor n/2\rfloor)$ is the largest value if $n$ is large enough.