One class of examples in a compact, connected Lie group $G$ of rank $r$ are the conjugacy classes of codimension $r$. Taking an element $a\in G$ whose centralizer is a maximal torus (a generic condition), the conjugacy class of $a$ is a submanifold $C_a\subset G$ of codimension $r$ whose normal plane at $a$ is the tangent at $a$ to $Z_a$, the centralizer of $a$ in $G$, which is a flat, totally geodesic submanifold. By the $G$-homogeneity of $C_a = G/Z_a\subset G$, the fact that this holds at $a$ implies that it holds at all point of $C_a$.
This gives an $r$-parameter family of mutually noncogruent examples.
There are other examples: When $G=\mathrm{SU}(n)$, and $a = \mathrm{diag}(\lambda_1,\lambda_2,\ldots,\lambda_n)$ is a diagonal element for which the $\lambda_i^2$ are all distinct, the submanifold
$$
M_a = \{\ g_1 a g_2\ |\ g_1,g_2\in\mathrm{SO}(n)\ \}\subset\mathrm{SU}(n)
$$
is homogeneous under the isometry group of $\mathrm{SU}(n)$ (endowed with its biïnvariant metric), and its tangent plane at $a$ is orthogonal to the diagonal maximal torus $T\subset\mathrm{SU}(n)$, so, by homogeneity, its tangent plane at any point is orthogonal to a flat, totally geodesic submanifold of $\mathrm{SU}(n)$. Hence it has an abelian normal bundle.
This gives another $r=n{-}1$ parameter family of mutually noncongruent examples of codimension $r$, distinct from the conjugacy classes.
Using the methods of exterior differential systems, one can show that, when $n=3$, these two families account for all of the codimension $2$ submanifolds of $\mathrm{SU}(3)$ with abelian normal bundle, in the sense that any connected submanifold $M^6\subset\mathrm{SU}(3)$ with abelian normal bundle is, up to ambient isometry, an open subset of one of the examples listed above. The argument that I have written out is a calculation using exterior differential systems and the moving frame, but, when I have time, I can sketch the proof, if there is interest.
Addendum 1: I had a flight with some time to look at the other two compact simple rank 2 groups. It turns out that all of the connected codimension two submanifolds with abelian normal bundle in $\mathrm{SO}(5)$ and $\mathrm{G}_2$ are (open subsets of) homogeneous compact ones as well. In each case, there is one additional $2$-parameter family of examples beyond the principal conjugacy classes.
Addendum 2: I also checked the rank $r=3$ case $G = \mathrm{SU}(4)$, and found that every connected codimension $3$ submanifold of $G$ with abelian normal bundle is an open subset of one of the two types of homogeneous examples listed above. Since any codimension $2$ submanifold of $\mathrm{SU}(4)$ that is foliated by codimension $3$ submanifolds of $\mathrm{SU}(4)$ with abelian normal bundle will have abelian normal bundle, it follows that there are many non-homogeneous codimension $2$ submanifolds of $\mathrm{SU}(4)$ with abelian normal bundle.
In fact, it now seems likely that, for any compact simple group $G$ of rank $r>1$, there are two $r$-parameter families of homogeneous codimension $r$ submanifolds with abelian normal bundle and every connected codimension $r$ submanifold with abelian normal bundle is, up to ambient isometry, an open subset of a homogeneous one. The two families are as follows: The first family is the family of principal conjugacy classes in $G$, and the second family is the set of principal orbits of $K\times K$ acting by left and right multiplication in $G$ where $G/K$ is a symmetric space of rank $r$. For example, when $G=\mathrm{SU}(n)$ $(n\ge3)$, $K=\mathrm{SO}(n)$; when $G = \mathrm{SO}(n)$ $(n\ge 5)$, $K = \mathrm{SO}(p)\times\mathrm{SO}(n{-}p)$ where $p = \bigl[\tfrac12 n\bigr]$; when $G=\mathrm{Sp}(n)$ $(n\ge3)$, $K=\mathrm{U}(n)$; when $G=\mathrm{G}_2$, $K=\mathrm{SO}(4)$; when $G=\mathrm{F}_4$, $K=\mathrm{Sp}(3)\mathrm{SU}(2)$; when $G=\mathrm{E}_6$, $K=\mathrm{Sp}(4)$; when $G=\mathrm{E}_7$, $K=\mathrm{SU}(8)$, and when $G=\mathrm{E}_8$, $K=\mathrm{SO}'(16)$.
Best Answer
It is known (say, Y. Tashiro, Complete Riemannian manifolds and some vector fields, Trans.Amer.Math.Soc. 117(1965) 251– 275; I am not sure that Tashiro is the first who proved it and there were many later papers which independently prove the same result later. The proof is pretty straightforward.) that the existence of such a function for a metric implies that the metric is a warped product metric, $$ g= dx_1^2 + \sigma(x_1) \sum_{i,j=2}^{n} h_{ij}(x_2,...,x_n)dx_i dx_j, \ \ \ (*) $$ and for such metrics a solution $f$ of your equation $\nabla\nabla f= \lambda(x) g$ is some function of $\sigma$.
The special case when the metric $g$ is flat is much more easy. In this case the function $\sigma(x_1)$ from $(*)$ is automatically $\textrm{const_1} \ (x_1+ \textrm{const}_2)^2$ (this fact is also classically known), so in the nontrivial case, when $\textrm{const_1}\ne 0$, after the affine reparameterization of $x_1$, the form $(*)$ is the cone form $$ g= dx_1^2 + x_1^2 \sum_{i,j=2}^{n} h_{ij}(x_2,...,x_n)dx_i dx_j, $$
This implies that in the flat case the function $\sigma(x_1)$, in some euclidean coordinate system $(y_1,...,y_n)$, is (after addition of linear terms and constant which change nothing) simply the function $$f= y_1^2 +...+ y^2_n.$$
For flat metrics of other signatures the answer is essentially the same, in this case $f= \varepsilon_1y_1^2 +...+ \varepsilon_ny^2_n,$ where $\varepsilon_i\in \{\pm 1\}$ are responcible for the signature