More specifically, $f$ harmonic s.t.:
$f: (0,1)^{2} \rightarrow \mathbb{R}$ with $f(x,y)=g(x)h(y)$ and $g,h \in C^{2}((0,1))$
Bear in mind that I am at calc III level, having just started an introductionary course for PDEs. This is our very first exercise. I have several ideas/starting points but things complicate so quickly that I have trouble estimating which one is promising and how to marry them.
First I deduced the obvious requirement of
$\forall x,y \in (0,1): g''(x)h(y)=-g(x)h''(y)$
Secondly, by mean value property, $f$ is of the given form iff
$\forall B_{r}((x,y)) \subseteq (0,1)^{2}: f(x,y)=g(x)h(y)=\frac{1}{\pi r^{2}}\int_{B_{r}((x,y))}g(x)h(y)d\mathcal{L}^{2}$
Where to go from here? Do I need to make specific computations and/or coordinate transformations with that integral? One idea I had was that if $f$ satisfies the above it is harmonic iff
$\frac{\partial^{2}}{\partial x^{2}}\int_{B_{r}((x,y))}g(x)h(y)d\mathcal{L}^{2}=-\frac{\partial^{2}}{\partial y^{2}}\int_{B_{r}((x,y))}g(x)h(y)d\mathcal{L}^{2}$
If I can verify whether all conditions to exchange order of differentiation and integration are met (or restrict $f,g,h$ such that they are), I could proceed with
$\int_{B_{r}((x,y))}\frac{\partial^{2}}{\partial x^{2}}g(x)h(y)d\mathcal{L}^{2}=-\int_{B_{r}((x,y))}\frac{\partial^{2}}{\partial y^{2}}g(x)h(y)d\mathcal{L}^{2} \iff \int_{B_{r}((x,y))}g''(x)h(y)d\mathcal{L}^{2}=-\int_{B_{r}((x,y))}g(x)h''(y)d\mathcal{L}^{2}$
In retrospect, these steps look like circular reasoning, I basically just integrated both sides of the equation I deduced at the start over arbitrary Balls s.t. $B_{r}((x,y)) \subseteq (0,1)^{2}$.
At this point I am a bit lost. Product rule? Change of coordinates/Fubini? I even have trouble knowing when I am done, because I have no idea how much you can narrow down the class of functions $g,h$ satisfying this.
Fix $y.$ Then your equation tells you that $h^{\prime\prime}(x) = c h(x),$ so on horizontal lines your function is a complex exponential in $x$ with coefficient depending on $y.$ The same is true on vertical lines, so the log of your function is linear in both $x$ and $y.$