In the Euclidean space (or Lorentz spacetime, if you are interested in relativity), there is one orthonormal coordinate system $\{x^\mu\}$ such that the distance squared is given by
$ds^2=(dx^1)^2+\cdots+(dx^n)^2$
I want to find a general coordinate transformation $y^\mu=\phi^\mu(x)$ such that the distane is
$ds^2=(dy^1)^2+\cdots+(dy^n)^2$
Being general, I mean $\partial\phi^\mu/\partial x^\nu$ are functions of $x^\mu$.
We know that we can find a particular coordinate transformation to make this happen, that is, rotation. Now, if we want someting more general, I guess the conditions are
$\frac{\partial\phi^\mu}{\partial x^\rho}\frac{\partial\phi^\nu}{\partial x^\sigma}g^{\rho\sigma}=g^{\mu\nu}$
and
$\nabla^2y^\mu=0$
The first shows that we can view $\frac{\partial\phi^\mu}{\partial x^\rho}$ as the $\rho$-th compoent of a vector labeled by $\mu$, and these vectors form an orthonormal basis. We can find infinitely many sets of such kind of orthnormal basis smoothly. But it does not seem to me they definitely correspond to a coordinate transformation probably because of the second condition. However, the second condition does not seem to impose any strong constraints on what kind of orthonormal basis we shall choose.
So I get stuck. Please help me out!
An obvious way to think about this is to say that since the metrics are the same, the geodesics are the same (viz. straight lines). Therefore the transformation has to map all straight lines to straight lines. Straight lines are of the form $$l^{\mu}(t) = a^{\mu}t+b^{\mu}$$ where $a^{\mu},b^{\mu}$ are constants, and it is easy to see that the only way to make this work is to take $\phi^{\mu} = R^{\mu}_{\nu}x^{\nu}+A^{\mu}$ with $R$ a constant matrix (i.e. $\phi$ is a linear map). Then your first condition gives that $R$ is a rotation (or a Lorentz transformation in the case of Minkowski space).