What's wrong with this tensor product approach to diagonalize a 4x4 matrix? I have a matrix that looks like:
\begin{align}
     H &= 
    \begin{bmatrix}
0 & 2\sin(2K_y) & 2\sin(2K_x) & 0\\
2\sin(2K_y) & 0 & 0 & \gamma\\
2\sin(2K_x) & 0 & 0 & 2\sin{2K_y}\\
0 & \gamma & 2\sin{2K_y} & 0\\
\end{bmatrix}
\end{align}
where $\gamma = 2\sin(2K_x + 2\pi/N)$, and $N = 2,3$. I write this in the tensor product form:
\begin{align}
H = \sigma^x \otimes \Bigg[(\frac{I_{2x2}+\sigma^z}{2})2\sin (2Kx) + (\frac{I_{2x2}-\sigma^z}{2}) \gamma\Bigg] + 2\sin(2K_y) I_{2x2} \otimes \sigma^x 
\end{align}
Let's say our solution looks like $\psi_x \otimes \psi_y$. And then suppose that we have a sub equation $\sigma^x \psi_x = l \psi_x$, then we can solve the eigenvalues of H by subbing this equation in and getting:
\begin{align}
\Bigg[l \Big[(\frac{I_{2x2}+\sigma^z}{2})2\sin (2Kx) + (\frac{I_{2x2}-\sigma^z}{2}) \gamma\Big] + 2\sin(2K_y) \sigma^x \Bigg] \psi_y = \lambda \psi_y
\end{align}
Of course, we have to make sure that the $\psi_x$ eigenequation has a solution, and it does for $l=\pm 1$.
We then solve the vector equation for $\psi_y$ to get the following eigenvalues:
\begin{align}
\lambda_{\pm} = \frac{l (2x + \gamma)}{2} \pm \frac{\sqrt{(2x + \gamma)^2 + 16y^2 - 8\gamma x}}{2}
\end{align}
Where I've used the short form $x = \sin (2K_x), y = \sin(2K_y)$.
Which gives me the following graph for N=2 (makes sense, as $\gamma \rightarrow -2\sin(2K_x)$ for N=2, thus the solution becomes $\lambda = \pm 2\sqrt{\sin^2 x + \sin^2 y}$):

And the following graph for N=3:

To check this, we can diagonalize the 4x4 matrix directly, without reducing it down to the tensor product form. It agrees with the tensor solution for N=2, but gives me a completely different graph for N=3! (It doesn't seem very apparent in this figure but even here we can see that they have different answers at the edges. And this has a depression that the tensor solution does not) 
I trust the straight diagonalization more than my tensor solution, but I can't place any obvious mistakes in my calculations. Is my tensor decomposition of H at fault? Is the agreement of the 2 graphs for the case of N=2 just a fluke?
 A: I wasn't entirely able to follow your approach, so here is another tensor product approach that might work for you. As you have noted, we have
$$
H = \sigma_x \otimes A + \alpha \cdot I\otimes \sigma_x
$$
for some matrix $A$ and constant $\alpha$. Denote
$$
U = \frac 1{\sqrt{2}}\pmatrix{1&1\\1&-1}.
$$
Note that $U$ is unitary and $U^\dagger \sigma_x U = \sigma_z$ (where $\dagger$ denotes a conjugate transpose). With that, we find that
$$
M := (U \otimes I)^\dagger H(U \otimes I) = (U^\dagger\sigma_xU) \otimes A + \alpha \cdot (U^\dagger U) \otimes \sigma_x
\\ = \sigma_z \otimes A + \alpha \cdot I \otimes \sigma_x.
$$
Notably, the matrix $M$ is block diagonal with
$$
M = \pmatrix{\alpha \,\sigma_x + A & 0\\0 & \alpha \ \sigma_x - A}.
$$
It now suffices to separately diagonalize these $2 \times 2$ matrices. In particular, suppose that $V^\dagger(\alpha \sigma_x + A)V$ and $W^\dagger (\alpha \sigma_x - A)W$ are diagonal. It follows that the matrix $P^\dagger M P$ is diagonal, with
$$
P = \pmatrix{V & 0\\0 & W}.
$$
Putting that all together, we conclude that
$$
\left[\pmatrix{V & 0\\0 & W} (U \otimes I)\right]^\dagger H \left[\pmatrix{V & 0\\0 & W} (U \otimes I)\right]
$$
is diagonal, which is to say that the columns of the matrix $\pmatrix{V & 0\\0 & W} (U \otimes I)$ are the eigenvalues of $H$.
