An effective/mechanical way to check if a function is separable By $f(x, y)$ being separable, I mean $f(x,y)$ can be expressed as $f(x,y)=g(x)h(y)$ for some functions $g$ and $h$.
By an effective procedure, I mean mechanical/mindless/automatic way to do this, one that can be articulated as a short algorithm, that conclusively determines one way or the other after some time whether $f$ is separable or not.
Would the answer be different if we considered only analytic/named functions, versus if we considered all functions (including only numerically defined, unnamed ones)?
 A: This problem is undecidable due to Rice's theorem.
A: It depends very much how $f$ is expressed. If it's a black box - input arguments and get values, no idea how or even if there's any sort of formula - then no, there isn't any way to find that it's separable. What there is is a way to have a very good chance of quickly finding when it's not separable, and thus being able to hypothesize that a function likely is separable if we go long enough.
First, I'd like to change the terminology; instead of saying "separable", we ask whether the tensor rank of a function is $1$. What's the tensor rank? It's the smallest possible $n$ such that we can write $f(x,y)=\sum_{i=1}^n g_i(x)h_i(y)$. For example, the tensor rank of $f(x,y)=0$ is $0$, the tensor rank of $f(x,y)=x^2\sin(y)$ is $1$, the tensor rank of $f(x,y)=(x+y)^2=x^2+2xy+y^2$ is $3$, and the tensor rank of $f(x,y)=\frac1{x+y+1}$ for $x,y\ge 0$ is $\infty$.
What good does this do? Well, we can also talk about the tensor rank of a matrix; the tensor rank of $A$ is the smallest $n$ such that we can write $A=\sum_{i=1}^n u_iv_i^T$, where the $u_i$ are column vectors and the $v_i^T$ are row vectors. If we build a matrix with $jk$ entry $f(x_j,y_k)$, the tensor rank of the function restricted to those points is the same as the tensor rank of the matrix - because we can exactly match a rank $1$ tensor $g_i(x)h_i(y)$ to a column by row product $G_i H_i^T$. The tensor rank of a general function (on a domain that's a Cartesian product) is equal to the supremum of the finite values we get from matrices of this form.
And then, there's a theorem: the tensor rank of a matrix is the same as the row rank, the column rank, or any other way we have to calculate it. There is just one rank for matrices, and this apparently new "tensor rank" is just another way to calculate rank.
So then, we have our method: Choose some random finite collection of values $x_j$ and $y_k$, find $f(x_j,y_k)$ for each $j$ and $k$, and look at whether the matrix $A$ of values $a_{jk}=f(x_j,y_k)$ has rank at most $1$. If it does, $f$ might have tensor rank $1$, and we can check again with more values. If it doesn't, $f$ has tensor rank greater than $1$ and we're done.
Why is it impossible to guarantee a decision? Consider $f(x,y)=\begin{cases} 2&(x,y)=(a,b)\\ 1&\text{otherwise}\end{cases}$ for $0\le x,y\le 1$ in $\mathbb{R}^2$. If our lattice of points doesn't include the one bad point $(a,b)$, all of the values will be $1$ and the matrix will have rank $1$. If it does include the bad point, the matrix will have rank $2$. The function has tensor rank $2$ - but our probability of finding it with any finite collection of test points is zero.
This is a worst-case scenario. Most of the time, when a function has tensor rank greater than $1$, we'll find that out quickly. Still, in order to positively prove that the tensor rank is $1$ in a case in which the domain isn't finite, we'll need to apply some other method using more knowledge of the function.
A: This is not meant to be rigorous, only suggestive of a blackbox test.
If g(x,y)=f(x)*h(y) then the ratio g(x,y1)/g(x,y2) does not depend on x. Pick values for y1 and y2 and make sure that neither is 0. Compare the ratio at different values of x -- also not zero -- and if they are the same then the function is (probably?) separable in x. Demonstration using SymPy:
>>> g = lambda x,y:cos(x+y)+cos(x-y)
>>> r = lambda x,y1,y2: (g(x,y1)/g(x,y2)).n()
>>> y1=2; y2=3
>>> r(1, y1,y2)
0.420353525886466
>>> r(2, y1,y2)
0.420353525886466

Both values are the same suggesting that g(x,y) is separable (and SymPy can identify that):
>>> simplify(g(x,y))
2*cos(x)*cos(y)

This is similar to the approach of testing to see if f'x/f is independent of y but doesn't compute the derivative. In both cases you will need to supply values for x and y since it will not -- until the separation has been identified -- be obvious that the equation is independent of y, e.g.:
>>> g(x, y).diff(x)/g(x, y)
(-sin(x - y) - sin(x + y))/(cos(x - y) + cos(x + y))
>>> _.is_constant(y)  # lots of internal expression wrangling...
True

See also notes here.
