We have the matrix eigenvalue problem
$$AX+XB^\text{T}=\lambda (CX+XD^\text{T})$$
Where $\lambda$ is the eigenvalue, $X$ is $m\times n$ and plays the role of the eigenvector, and $A$, $B$, $C$, and $D$, are square matrices of appropiate size, and for my purposes, $C$ and $D$ are diagonal. Rearranging we get to
$$(A-\lambda C)X=-X(B-\lambda D)^\text{T}$$
Now, trying to solve this through separation of variables, choose $X=uv^\text{T}$
$$(A-\lambda C)uv^\text{T}=-uv^\text{T}(B-\lambda D)^\text{T}=muv^\text{T}$$
where $m$ is the separation constant. Now we separate
$$(A-\lambda C)u=mu$$ $$(B-\lambda D)v=-mv$$
My question is, how can I solve this system of eigenproblems? Is there a numerical algorithm for this?
Background: all this came from the matrix analogy of solving a PDE of the form
$$ \frac{1}{f(x)+g(y)}(L_x u+ L_y u)=\lambda u $$
Where $L_x$ and $L_y$ are univariate linear differential operators.