Let $\Omega$ be a bounded domain in $\mathbb{R}^3$ and $f_1,f_2 \in C^2(\bar{\Omega})$. Suppose
$\int_{\Omega}(f_2-f_1)\varphi \, dx=0$ and $\int_{\Omega}(f_2 \Delta^{-1} f_2- f_1 \Delta^{-1} f_1)\varphi \, dx =0$
for all harmonic functions $\varphi$, where $\Delta^{-1}f=\int_{\Omega}f(y)\Phi(x-y) \, dy$ and $\Phi$ is the fundamental solution of the Laplacian. Is $f_1=f_2$?
Assume without loss of generality that $0\in \Omega$.
Let $f_1, f_2$ be radial, with support contained in $\Omega$.
By the mean value property of harmonic functions, you have
$$ \int f_2 \varphi ~\mathrm{d}x = 4\pi \varphi(0) \int_0^{\infty} r^2 f_2(r) ~\mathrm{d}r = \varphi(0) \int f_2(r) ~\mathrm{d}x. $$
Now let $f_1, f_2$ be arbitrary non-trivial radial functions such that $\int f_1 ~\mathrm{d}x = \int f_2 ~\mathrm{d}x= 0$.
As $\Delta^{-1} f_1$ and $\Delta^{-1}f_2$ are also radial, we have that $$\int f_1 \Delta^{-1} f_1 \varphi ~\mathrm{d}x = \varphi(0) \int f_1 \Delta^{-1} f_1 ~\mathrm{d}x$$ whenever $\varphi$ is harmonic. So then by replacing $f_1 \mapsto \lambda_1 f_1$ and $f_2 \mapsto \lambda_2 f_2$ we can certainly arrange for
$$ \int (\lambda_1)^2 f_1 \Delta^{-1} f_1 - (\lambda_2)^2 f_2 \Delta^{-2} f_2 ~\mathrm{d}x = 0 $$
and so as long as our seeds $f_1$ and $f_2$ have, for example, unequal supports, we arrive at a counterexample.
In the comments below the OP asked for examples where $f_1$ and $f_2$ are signed. Hopefully the computations below are (more or less) correct:
Roughly speaking, let $1 > r > \rho$ be fixed numbers. The idea is to let $f_1$ be the sum of the surface measure of the sphere of radius 1 with $a_\rho$ times the surface measure of the sphere of radius $\rho$, while $f_2$ is $a_r$ times the surface measure of the sphere of radius $r$. Then the first condition requires
$$ 1 + a_\rho \rho^2 = a_r r^2 $$
For the surface measure of a sphere of radius $\rho$, $\Delta^{-1}$ can be computed to be the function equaling $C |x|^{-1} \rho^2$ for some universal constant $C$. So you would have that the condition on $f_1 \Delta^{-1} f_1$ to give rise to the condition
$$ 1 + a_\rho^2 \rho + 2 a_\rho \rho^2 = a_r^2 r $$
$$ \implies 1 + a_\rho^2 \rho + 2 a_\rho \rho^2 = r^{-3} (1 + 2 a_\rho \rho^2 + a_\rho^2 \rho^4) $$
$$ \implies 1 - r^3 + 2 a_\rho \rho^2 (1 - r^3) + a_\rho^2 \rho (\rho^3 - r^3) = 0$$
Since $1 > r > \rho$ this last equation has a positive root $a_\rho$, which then gives a positive solution $a_r$.
To move from measures to functions: instead of using surface measures use $C^\infty_c$ functions highly concentrated near the corresponding surfaces. The computations would be slightly less explicit but the relationships between the coefficients $a_r, a_\rho$ would still hold up to a (small) perturbation, and you still can find a counterexample.
