Second, one can assume, by adding a constant to $u$, that $u(0,0)=0$, so I will assume this normalization made henceforth. Then, the obvious integral inequality arising from $|\nabla u| = |xy|^n$ and the Cauchy-Schwartz inequality would imply that $$ |u(x,y)| \le \frac{|xy|^n\sqrt{x^2+y^2}}{(2n{+}1)}. $$ In particular, $u$ would vanish to order $2n{+}1$ at $(0,0)$ and would satisfy $u(x,0) = u(0,y) = 0$. It follows from this that $u$ could not be of differentiability class $C^{2n+2}$, since, if it were, the limit function $$ p(x,y) = \lim_{r\to0} \frac{u(rx,ry)}{r^{2n+1}} $$ would exist and be a polynomial homogeneous of degree $2n{+}1$ that satisfied $|\nabla p|^2 = (xy)^{2n}$, and it is easy to show that there is no such polynomial. However, as will be seen, when $n\ge3$, there exists a $u\in C^{n-1}(\mathbb{R}^2)$ satisfying $|\nabla u|^2 = (xy)^{2n}$ and the homogeneity condition $u(rx,ry) = r^{2n+1}\,u(x,y)$$u(rx,ry) = |r|^{2n+1}\,u(x,y)$ for all $r$. This $u$ is real-analytic away from the lines $x\pm y = 0$ but fails to be $C^n$ on these two lines.
Now, computation shows that the Gauss curvature of $g$ is $K = n(x^2{+}y^2)/(xy)^{2n+2}>0$, which suggests that nearly all of the geodesics of $g$ will avoid going to the singular boundary where $xy=0$, and, indeed, this turns out to be the case (see below). To parametrize the geodesics, it turns out to be convenient to use a parameter $t$ other than arc length. A curve $\bigl(x(t),y(t)\bigr)$ in the first quadrant parametrizes a $g$-geodesic when there is a function $\phi(t)$ satisfying the ODE system $$ \dot x = xy\cos\phi, \quad \dot y = xy\sin\phi, \quad \dot\phi = n\,(x\cos\phi-y\sin\phi),\tag1 $$ and every $g$-geodesic in the first quadrant has such a parametrization, unique up to replacing $t$ by $t+t_0$ for some constant $t_0$. In this case, arclength $s(t)$ along the geodesic satisfies $\dot s = (xy)^n$$\dot s = (xy)^{n+1}$. [The advantage of writing the geodesic equations this way is that they extend smoothly across the singular locus $xy=0$.] Note that these equations are invariant under the homothetical scaling $(t,x,y,\phi)\to(t/r,rx,ry,\phi)$. Because of the scaling symmetry of the equations, one can extract a 2D phase portrait that makes clear the behavior of the geodesics as follows: Let $x+iy = \mathrm{e}^{u+iv}$. Then the above equations become (after a change of independent variable) $$ u' = \cos(\phi{-}v)\,\cos v\sin v,\qquad v' = \sin(\phi{-}v)\,\cos v\sin v,\qquad \phi' = n\,\cos(\phi{+}v).\tag2 $$
One can now draw the $v\phi$-phase portrait, concentrating on the rangestrip $0\le v\le \pi/2$, which represents the first quadrant in the $xy$-plane, and bearing in mind that thethese equations are invariant under the translationinvolution $(v,\phi)\to (v{+}\pi,\phi+\pi)$$(v,\phi)\to(\tfrac12\pi{-}v,\tfrac12\pi{-}\phi)$ and reverse under $(v,\phi)\to (v,\phi+\pi)$$(v,\phi)\to (v,\phi{+}\pi)$. There is
There are a sink at $S_- = (v,\phi)=(\pi/4,\pi/4)$, a source at $S_+ = (v,\phi)=(\pi/4,-3\pi/4)$, and saddles at $S_1 = (v,\phi)=(0,\pi/2)$, $S_2 = (v,\phi)=(\pi/2,0)$, $S_3=(v,\phi)=(0,-\pi/2)$, and $S_4=(\pi/2,-\pi)$. There In addition to the 'trivial' separatrices that make up the boundary lines $v=0$ and $v=\tfrac12\pi$, there are 4four 'non-trivial' separatrices: $L_1$ leaving $S_1$ and going to $S_-$, $L_2$ leaving $S_2$ and going to $S_-$, $L_3$ leaving $S_+$ and going to $S_3$, and $L_4$ leaving $S_+$ and going to $S_4$.