Consider the multivariate Stein equation $Δf - x\cdot \nabla f = h(x) - E h(Z)$, where $Z$ is a standard $d$-dimensional Gaussian random vector, and let $f\_h$ be the solution given by Barbour's generator approach. We prove that, when $h$ is $α$-Hölder ($0<α\leq1$), all derivatives of order $2$ of $f\_h$ are $α$-Hölder {\it up to a $\log$ factor}; in particular they are $β$-Hölder for all $β\in (0, α)$, hereby improving existing regularity results on the solution of the multivariate Gaussian Stein equation. For $α=1$, the regularity we obtain is optimal, as shown by an example given by Raič \cite{raivc2004multivariate}. As an application, we prove a near-optimal Berry-Esseen bound of the order $\log n/\sqrt n$ in the classical multivariate CLT in $1$-Wasserstein distance, as long as the underlying random variables have finite moment of order $3$. When only a finite moment of order $2+δ$ is assumed ($0<δ<1$), we obtain the optimal rate in $\mathcal O(n^{-\fracδ{2}})$. All constants are explicit and their dependence on the dimension $d$ is studied when $d$ is large.