Your solution for $2$ dimensions works for any dimension. You simply solve for co-ordinates $x^1, \dots, x^n$ such that the metric tensor $g_{ij}dx^idx^j$ satisfies $$ g_{11} = g_{22},\ g_{12} = 0 $$ and set $f = x^1$, $g = x^2$. This is effectively parameterized isothermal co-ordinates, so your title is justified. You should be able to get solutions near a point by fixing a foliation of $2$-dimensional surfaces and solving for isothermal co-ordinates along each surface in such a way that the solution depends smoothly on the other $n-2$ co-ordinates.