As in my previous solution in the 3-dimensional case (discussed here), we can set $f=-(p/h)x$ for a function $x$ that satisfies
$$
\Delta x + p\,x^2 = |\nabla x|^2 - x^3 = 0.\tag 1
$$
Conversely, if $x$ satisfies this system for a metric $g$, then $f = -(p/h)x$ will satisfy the original equations.
Again, the same argument as before shows that the $1$-form $\omega_1 = x^{-3/2}\, \mathrm{d}x$ is a $1$-form of $g$-norm equal to 1, so one can write $g$ locally in the form
$$
g = {\omega_1}^2 + {\omega_2}^2 + \cdots + {\omega_n}^2.
$$
Again, fixing the orientation so that $\omega_1\wedge\cdots\wedge\omega_n$ is the oriented volume form, we have $\ast_g\omega_1 = \omega_2\wedge\cdots\wedge\omega_n$, and the definition of the Laplacian gives us $\mathrm{d}(\ast_g\mathrm{d}x) = -\Delta x\,\,{\ast_g}1$. Just as before, we now compute, using $\mathrm{d}x = x^{3/2}\,\omega_1$ that
$$
\mathrm{d}(x^{3/2-p}\,\omega_2\wedge\cdots\wedge\omega_n) = 0,
$$
and hence, locally, one can write
$$
x^{3/2-p}\,\omega_2\wedge\cdots\wedge\omega_n = \mathrm{d}y^2\wedge\cdots\wedge\mathrm{d}y^n\tag2
$$
for some functions $y^2,\ldots,y^n$ that, with $x$, gives a local coordinate system. Moreover, (2) implies that
$$
{\omega_2}^2 + \cdots + {\omega_n}^2 = x^{(2p-3)/(n-1)}\,g_{ij}(x,y)\,\mathrm{d}y^i\mathrm{d}y^j,
$$
for some function $g_{ij}$ satisfying $\det(g_{ij}) = 1$.
Conversely, given functions $g_{ij}$ on a domain in $\mathbb{R}^n$ so that $(g_{ij})>0$ and $\det(g_{ij}) = 1$, we have that the metric
$$
g = x^{-3}\,(\mathrm{d}x)^2 + x^{(2p-3)/(n-1)}\,\bigl(g_{ij}(x,y)\,\mathrm{d}y^i\mathrm{d}y^j\bigr)
$$
and the function $x$ satisfy the equation (1).
Finally, the equation $S =-p(p+1) f$ is a second order PDE for the functions $g_{ij}$, and the existence of solutions should be fairly straightforward as long as $n>2$. It's not clear that it's worth writing out the details.
