Solutions to PDE $\langle\nabla\psi,\nabla\psi+\mathbf{f} \rangle= 0$ Consider the vector field
$$
\mathbf{f}(x_1,x_2)=\begin{bmatrix}f_1(x_1,x_2)\\ f_2(x_1,x_2)\end{bmatrix}=\begin{bmatrix}\frac{1}{1+x_2}+x_1\\ \frac{1}{1+x_1}+x_2\end{bmatrix}.
$$
I'm interested in finding a scalar function $\psi(x_1,x_2)$ such that:
$$\tag{$\ast$}\label{ast}
\langle\nabla\psi,\nabla\psi+\mathbf{f} \rangle= \frac{\partial\psi(x_1,x_2)}{\partial x_1} \left(f_1(x_1,x_2)+\frac{\partial\psi(x_1,x_2)}{\partial x_1}\right)+\frac{\partial\psi(x_1,x_2)}{\partial x_2} \left(f_2(x_1,x_2)+\frac{\partial\psi(x_1,x_2)}{\partial x_2}\right)=0
$$

My question. Does there exist a solution $\psi(x_1,x_2)$ to \eqref{ast}? If so, how to compute it?

I'm aware that my question could be a trivial or naive one, but since I'm new to this kind of problems I would really appreciate any help/comment/suggestion. Thank you.
 A: On the one hand, the nullity of the scalar product
$$
\nabla\psi\cdot (f+\nabla\psi) = \psi_{,1} (f_1+ \psi_{,1} )
+ \psi_{,2} (f_2+ \psi_{,2} ) = 0
$$
leads to the orthogonality of $\nabla\psi$ and $f+\nabla\psi$. Thus, $f+\nabla\psi$ is proportional to the vector $\nabla\psi^\perp = (\psi_{,1}, -\psi_{,2})^\top$. There exists $\alpha(x_1,x_2)$ such that $f+\nabla\psi = \alpha \nabla\psi^\perp$, or equivalently
$$
(1-\alpha)\psi_{,1} = -f_1
\qquad\text{and}\qquad
(1+\alpha)\psi_{,2} = -f_2 .
$$
By eliminating successively $\alpha$ and unity of the previous system, we get to
\begin{aligned}
2\psi_{,1}\psi_{,2} &= -f_2\psi_{,1} - f_1\psi_{,2}
\\
2\alpha\psi_{,1}\psi_{,2} &= -f_2\psi_{,1} + f_1\psi_{,2} .
\end{aligned}
On the other hand, the nullity of the scalar product rewrites as
$$
(\psi_{,1} + \psi_{,2})^2 - 2\psi_{,1}\psi_{,2} + f_1\psi_{,1} + f_2\psi_{,2} = 0 ,
$$
so that previous identities lead to
$$
(\psi_{,1} + \psi_{,2})^2 + (f_1+f_2)(\psi_{,1}  + \psi_{,2}) = 0 .
$$
Thus, the solutions satisfy $\psi_{,1} + \psi_{,2} = -(f_1+f_2)$ if they have nonzero divergence -- otherwise, they satisfy $\psi_{,1} + \psi_{,2} = 0$, which is a particular case of the previous non-homogeneous equation. Solutions to these linear first-order PDEs can be obtained by using the method of characteristics for the Lagrange-Charpit system
$\text d x_1 = \text d x_2 = {\text d \psi}/{r}$ with $r = -(f_1+f_2)$ (see related posts on this site). Solutions are of the form
$$
\psi(x_1, x_2) = \int^{x_1} r(\xi,\xi + x_2-x_1)\, \text d \xi + F(x_2-x_1) ,
$$
where $F$ is an arbitrary function. The latter is determined by injecting the previous expression in the orthogonality condition $\nabla\psi\cdot (f+\nabla\psi) = 0$.
If $f$ is constant, then linear solutions of the form
\begin{aligned}
\psi(x_1,x_2) &= -f_1x_1 - f_2x_2 + C \\
\text{or}\qquad
\psi(x_1,x_2) &= -\tfrac12(f_1+f_2)(x_1+x_2)+ C 
\end{aligned}
are obtained.
In particular, one notes that the solutions are constant if $f \equiv 0$.
In theory, such solutions can be obtained similarly in the case
$$
f_1(x_1,x_2) = x_1 + (1+x_2)^{-1}, \qquad
f_2(x_1,x_2) = x_2 + (1+x_1)^{-1} ,
$$
but computations are more involved.
