This result is well-known in the theory of harmonic morphisms, about which, there is an extensive literature. It is a quite general fact (not depending on the conformally flat case of Euclidean space), implying that, when $n>2$, any conformal map $f:(M^n,g)\to (N^n,h)$ between (connected) Riemannian manifolds of dimension $n$ that is harmonic is necessarily a homothety i.e., $f^*(h) = r^2 g$ for some constant $r\ge0$.
You can find proofs in introductory papers and books on harmonic morphisms. I'm traveling, so I don't have time to check the best references, but, for example, try B. Fuglede, Harmonic morphisms between Riemannian manifolds, Ann. Inst. Fourier 28 (1978), 107-144, which probably has a proof written out there.
You can also find a proof in the introductory sections of my paper, Harmonic morphisms with fibers of dimension one (Communications in Analysis and Geometry, Volume 8 (2000), 219–265, but you'll have to sort out the relevant bits of the discussion in order to trace through a proof.