Let us reformulate OP's question as follows:
Give a proof that a local coordinate transformation $x^{\mu} \to y^{\rho}=y^{\rho}(x)$ between two local coordinate systems (on a 3+1 dimensional Lorentzian manifold) must be affine if the metric $g_{\mu\nu}$ in both coordinate systems happen to be on constant flat Minkowski form $\eta_{\mu\nu}$.
Here we will present a proof that works both with Minkowski and Euclidean signature; in fact for any signature and for any finite non-zero number of dimensions, as long as the metric $g_{\mu\nu}$ is invertible.
1) Let us first recall the transformation property of the inverse metric $g^{\mu\nu}$, which is a contravariant $(2,0)$ symmetric tensor,
$$ \frac{\partial y^{\rho}}{\partial x^{\mu}} g^{\mu\nu}_{(x)}\frac{\partial y^{\sigma}}{\partial x^{\nu}}~=~g^{\rho\sigma}_{(y)}, $$
where $x^{\mu} \to y^{\rho}=y^{\rho}(x)$ is a local coordinate transformation. Recall that the metric $g_{\mu\nu}=\eta_{\mu\nu}$ is the flat constant metric in both coordinate systems. So we can write
$$ \frac{\partial y^{\rho}}{\partial x^{\mu}} \eta^{\mu\nu}\frac{\partial y^{\sigma}}{\partial x^{\nu}}~=~\eta^{\rho\sigma}. \qquad (1) $$
2) Let us assume that the local coordinate transformation is real analytic
$$y^{\rho} ~=~ a^{(0)\rho} + a^{(1)\rho}_{\mu} x^{\mu} + \frac{1}{2} a^{(2)\rho}_{\mu\nu}x^{\mu}x^{\nu} + \frac{1}{3!} a^{(3)\rho}_{\mu\nu\lambda}x^{\mu} x^{\nu} x^{\lambda} + \ldots. $$
By possibly performing an appropriate translation we will from now on assume without loss of generality that the constant shift $ a^{(0)\rho} =0 $ is zero.
3) To the zeroth order in $x$, the equation $(1)$ reads
$$ a^{(1)\rho}_{\mu} \eta^{\mu\nu}a^{(1)\sigma}_{\nu}~=~\eta^{\rho\sigma}, $$
which not surprisingly says that the matrix $a^{(1)\rho}_{\mu}$ is a Lorentz (or an orthogonal) matrix, respectively. By possibly performing an appropriate "rotation", we will from now on assume without loss of generality that the constant matrix
$$ a^{(1)\rho}_{\mu}~=~\delta^{\rho}_{\mu} $$
is the unit matrix.
4) In the following, it will be convenient to lower the index of the $y^{\sigma}$ coordinate as
$$y_{\rho}~:=~\eta_{\rho\sigma}y^{\sigma}.$$
Then the local coordinate transformation becomes
$$y_{\rho} ~=~ \eta_{\rho\mu} x^{\mu} + \frac{1}{2} a^{(2)}_{\rho,\mu\nu}x^{\mu}x^{\nu} + \frac{1}{3!} a^{(3)}_{\rho,\mu\nu\lambda}x^{\mu} x^{\nu} x^{\lambda}+ \ldots$$
$$+\frac{1}{n!} a^{(n)}_{\rho,\mu_1\ldots\mu_n}x^{\mu_1} \cdots x^{\mu_n}+ \ldots. $$
5) To the first order in $x$, the equation $(1)$ reads
$$ a^{(2)}_{\rho,\sigma\mu}+a^{(2)}_{\sigma,\rho\mu}~=~0.$$
That is, $a^{(2)}_{\rho,\mu\nu}$ is symmetric in $\mu\leftrightarrow \nu$, but antisymmetric in $\rho\leftrightarrow \mu$. It is not hard to see (by applying the symmetry and the antisymmetry property in alternating order three times each), that the second order coefficients $a^{(2)}_{\rho,\mu\nu}=0$ must vanish.
6) To the second order in $x$, the equation $(1)$ reads
$$ a^{(3)}_{\rho,\sigma\mu\nu}+a^{(3)}_{\sigma,\rho\mu\nu}~=~0.$$
That is, $a^{(3)}_{\rho,\mu\nu\lambda}$ is symmetric in $\mu\leftrightarrow \nu\leftrightarrow \lambda $, but antisymmetric in $\rho\leftrightarrow \mu$.
For fixed $\lambda$, we can again reach the conclusion $a^{(3)}_{\rho,\mu\nu\lambda}=0$.
7) Similarly, we conclude inductively that the higher order coefficients $a^{(n)}_{\rho,\mu_1\ldots\mu_n}=0$ must vanish as well. So $y^{\mu}= x^{\mu}$. Q.E.D.