I have two matrix equations:
$$\begin{cases} XA_1=B_1Y \\ XA_2=B_2Y \end{cases}$$
where $A_1$, $B_1$, $A_2$, $B_2$ are known invertible square $n\times n$ matrices (and not the degenerate case where $A_1=kA_2$ annd $B_1=kB_2$ where $k$ is some constant), and $X$ and $Y$ are two square $n\times n$ matrices that I would like to solve for.
Is there any way to solve this "system of matrix equations," without explicitly writing out all $2n^2$ scalar equations? Is there any way to solve such equations in MATLAB or Mathematica?
In my specific use case, there can be more than two pairs of $A$, $B$ matrices if that helps, but I feel like that's just over constraining the problem.
The equations give $X=B_1YA_1^{-1}$ with $B_1YA_1^{-1}=B_2YA_2^{-1}$. So it is enough to solve the last equation. It is of the type $AYB-CYD=E$, which has been studied in many papers, e.g., here, or here.