Let $\Omega \subset \mathbb R^p$ be a convex, bounded domain with a smooth boundary. Let $a_{ij} : \Omega \to \mathbb R_+$ be a non-negative smooth function for $i, j \in \{1, 2\}.$ I am interested in the following system of PDE in terms of $f_1, f_2: \Omega \to \mathbb R$: $$ \begin{cases} \text{div} (a_{11} \nabla f_1 + a_{21} \nabla f_2) = 0 \text{ on } \Omega\\ \text{div} (a_{12} \nabla f_1 + a_{22} \nabla f_2) = 0 \text{ on } \Omega\\ \langle a_{11} \nabla f_1 + a_{21} \nabla f_2, n \rangle = 0 \text{ on } \partial \Omega\\ \langle a_{12} \nabla f_1 + a_{22} \nabla f_2, n \rangle = 0 \text{ on } \partial \Omega.\\ \end{cases} $$
If matrix $(a_{ij})$ is positive definite (even if it is not symmetric), then it is easy to verify that $f_1 = f_2 = \text{cnst.}$ is the unique solution. I am wondering if this positive definiteness condition can be relaxed while maintaining the uniqueness. In particular, is the invertibility of the $(a_{ij})$ sufficient for ruling out non-trivial solutions?
If $a_{ij}$ is constant, the invertibility implies $\Delta f_i = 0$ and $\langle \nabla f_i, n \rangle = 0,$ so the solution is unique. I conjecture that when $a_{ij}$ varies, the uniqueness does not hold in general, but I cannot find a counterexample.
