En primer lugar, la EDP puede hacerse homogénea dejando que $f=\exp{g(x,y)}$ donde $g(x,y)$ debe divergir a $-\infty$ como $\|(x,y)\|\to\infty$ de forma que se cumpla la condición de $f$ . Reescritura de la EDP, \begin{equation} \frac{\partial g}{\partial x}\frac{\partial g}{\partial y}+2\frac{\partial^2 g}{\partial x\, \partial y}=0. \end{equation}
En segundo lugar, dado que la EDP es simétrica con respecto al intercambio entre $x$ y $y$ su solución $f(x,y)$ también debe ser simétrica, al igual que $g(x,y)$ . Esta simetría implica que las derivadas parciales de $g(x,y)$ con respecto a $x$ y $y$ son iguales, \begin{equation} \frac{\partial g(x,y)}{\partial x}=\frac{\partial g(y,x)}{\partial y}=\frac{\partial g(x,y)}{\partial y}. \end{equation}
Dado que las derivadas parciales son las mismas, $g(x,y)$ puede dejarse igual a una función de una sola variable $h(r)$ donde $\frac{\partial r}{\partial x}=\frac{\partial r}{\partial y}=1$ es decir $r=x+y$ . Por tanto, la solución de la EDP está en la solución de \begin{equation} \left(\frac{\partial h}{\partial r}\right)^2 + 2 \frac{\partial^2 h}{\partial r^2}=0. \end{equation}
Desgraciadamente, no creo que exista una solución a esta ecuación que satisfaga $h(r)\to -\infty$ como $r\to\infty$ . Si se suprime esta condición, la solución es de la forma \begin{equation} h(r)=2\log({|r+C_1|})+C_2, \end{equation} que al volver a sustituirlo da una solución similar a la que el OP ha calculado, \begin{equation} f(x,y)=(A(x)+A(y))^2. \end{equation}
Espero que este trabajo no concluyente ayude.