Criterion of function separability in which function is plugged in differential equation. Question.
What is procedure of determining whether function is separable?
E.g., there are $g$, $h$ such that $$f(x,y)=g(x)h(y) \tag{1}$$
Is the following procedure valid? If so, why?

Relevance.
I was trying to understand why wave-function in Schrodinger equation describing motion of particle in two-dimensional box is separable.
$$ {\frac{-\hbar^2}{2m}}\left(\frac{\partial^2 f(x,y)}{\partial x^2} +\frac{\partial^2 f(x,y)}{\partial y^2} \right)=E f(x,y)~~~(A1)$$
I have met similar issues trying to prove separability of the function describing motion of particle rotating in spherical shell and separability of the function describing electron motion around hydrogen atom.
Method validity which I want to prove.
As far I understood the approach in the textbook I was reading, to prove that function is separable in is sufficient to take the differential equation and plug inside it product of two functions. If it is possible rearrange the result of this substitution into a sum of two parts (that is equal to constant), where each of them contains only one variable, than the function is separable. I have trouble understanding why this approach is valid.
Here is example from Atkins Physical Chemistry textbook.
So by substituting (1) into (A1) we obtain this.
$$ {\frac{-\hbar^2}{2m}}\left(\frac{\partial^2 [g(x)h(y)]}{\partial x^2} +\frac{\partial^2 [g(x)h(y)]}{\partial y^2} \right)=E [g(x)h(y)]~~~(A2)$$
$$ {\frac{-\hbar^2}{2m}}\left(\frac{d^2 [g(x)]}{d x^2} +\frac{d^2 [h(y)]}{d y^2} \right)=E [g(x)h(y)]~~~(A3)$$
$$ {\frac{-\hbar^2}{2m}}\left(\frac{d^2 [g(x)]}{d x^2} \right)=E_x [g(x)]~~~(A4)$$
$$ {\frac{-\hbar^2}{2m}}\left(\frac{d^2 [h(y)]}{d y^2} \right)=E_y [h(y)]~~~(A5)$$
So it is concluded that f(x,y) is indeed separable. By similar argument it can be proven that function describing motion particle in spherical shell is separable and also wavefunction describing electron motion around hydrogen atom is separable.
My attempt to establish validity of this method.
I know that any differential equation with two variables can be expressed as
$$m(x,y,z(z,y), f'x,f'y, f''xy... )=0 ~~~(2)$$
I have tried expressing the fact that function is separable without quantification over functions, because I don't know how to work with them. So I arrived at
$$\forall x \forall y [f(x,y) f(y,x)=f(x,x) f(y,y)] ~~~(3) $$
I obtained this equation by making substitutions x=y in the (1) equation and combining results with it again. This equation follows from the separability of functions.
I also obtained this equation.
$$\forall a\forall b\forall x \forall y [f(x,y) f(b,a)=f(x,a) f(b,y)] ~~~(4) $$
It can be obtained by making substitution x=a and y=b in (1) equation, combining both equations and using (1) equation again. Obtained equation should be equivalent with the fact that f(x,y) is separable.
I suspect that differential equation where f(x,y) is substituted with g(x)h(y) can be expressed as sum of equations each of which contains only one variable can be expressed as follows.
$$ m(x,y,z(x,y), f'x,f'y, f''xy... )=0\iff m1(x,y,z1(x), f'x, f''x... )+m1(x,y,z2(y), f'y, f''y... )+h=0~~~(5)$$
So I guess I need to prove that (5) implies (1).
In that textbook I was reading this was presented as something as self-evident, so I am grasping that answer should be simple.
Maybe it has something to do with fact that these differential equations are linear operators.
My second attempt to solve this question.
I tried to plug in some easy to handle inseparable function like this. $$f(x,y)=x^2+y^2 \tag{B1}$$
$$\frac{\partial^2 (x^2+y^2)}{\partial x^2}=2 \tag{B2}$$
$$\frac{\partial^2 (x^2+y^2)}{\partial y^2}=2 \tag{B3}$$
So if we plugg previous equations in (A1) we get this.
$$ {\frac{-\hbar^2}{2m}}\left(2+2 \right)=E(x^2+y^2) \tag{B4}$$
Which can be simplified to this.
$$x^2+y^2=-\frac{2\hbar}{Em} \tag{B5}$$
Which is valid equation, so it must be a solution of (A1), so if my reasoning is correct, this means that inseparable functions can be solution of (A1).
My third attempt to solve this question.
So in last attempt I didn't consider boundary conditions. Now lets take this function.
$$f(x,y)=x^2+y^2-k^2\tag{C1}$$, where $k2$ is some constant.
By similar procedure as in last attempt after plugging in this function inside (A1) we obtain this.
$$x^2+y^2=\frac{-2\hbar^2}{m E}+k^2 \tag{C2}$$
Now lets take that boundary conditions are
$$x^2+y^2>\gamma^2 \implies f(x,y)=0\tag{C3}$$, where $\gamma$ is some numeric constant.
But our function when goes out some circle of radius $\gamma$ is negative. So we redefine it as peace-wise function
$${ \begin{cases}
f(x,y)=x^2+y^2-k^2, x^2+y^2<=\gamma \\
f(x,y)=0 ,x^2+y^2>\gamma 
\end{cases} } \tag{C4}$$
Now for case when $x^2+y^2<\gamma^2$ solution of (A1) is the same.
For case when $x^2+y^2>\gamma^2$ soluition will be this.
$$ {\frac{-\hbar^2}{2m}}\left(0+0 \right)=E(0) \tag{C5}$$
$$0=0$$
So this equation satisfies (A1).
So if my reasoning correct if particle is in 2D-box with circular walls, than solutions can be inseparable.
P.S. I am a chemist. Mathematicians and physicians often find some things that other people don't understand as self evident. They are self evident for them probably because they work with them frequently. This why proofs are needed. So I need list of statements that starts from things I believe and than by known inference rules conclusion is derived. Each step must be executed formally. It is good if proof comes from book or publication.
 A: Your initial question is sort of different from the context you provided in Relevance. You see, you asked initially:

What is the procedure for determining whether a function, say $f(x,y)$, is separable?

But then you said:

I was trying to understand why the wave-function in the Schrödinger equation describing the motion of a particle in a two-dimensional box is separbale.

These are sort of two different questions, though they are related. I'm pointing this out because I need to organize my answer in a way that is clear to you & addresses both of these questions.
Okay, so there is no general procedure for determining whether a given function is separable in the sense that you described. It depends on the function and usually, you'd have to do some kind of extra (algebraic) manipulation to get it into separated form. For instance:
$$f(x,y) = x\sin(y)$$
is clearly separable while:
$$f(x,y) = \exp(-\alpha x) \cos(y) + 2\exp(-\alpha x) \sin(y)$$
is separable but you need to factorize in order to write it in a separable form. Once again, there isn't a general procedure for doing this with an arbitrary function.
Now, as for your question regarding wave-functions, the thing to understand is that people usually tend to assume that we are looking for separable solutions. So, we assume that a solution of the form:
$$f(x,y) = g(x)h(y)$$
satisfies the Schrödinger equation, do some formal trickery to find an explicit form for a (possible) solution and then proceed to check if that candidate solution actually works. Once all of that is done, we then proceed to try and prove that the solution is unique.
The way in which we would prove that a solution is unique depends on the problem. There are definitely ways to do that and it would be much better for me to just direct you to a book on PDEs which will surely contain an exposition of these methods.
I'm hoping that the above is sufficient to address the issues that come up later in your post. I just want to make a comment on your example from Atkins's book. So, given what I said above, you assume that we are looking for separable solutions. That is:
$$f(x,y) = g(x)h(y)$$
Then, it follows that:
$$\partial_x^2 f(x,y) = g''(x) h(y) \ , \ \partial_y^2 f(x,y) = g(x)h''(y)$$
So, in fact, the Schrödinger equation can be written as:
$$-\frac{\hslash^2}{2m} \left(g''(x)h(y) + g(x)h''(y) \right) = E(g(x)h(y))$$
