I don't have a completely general method for doing what you ask. However, there are a few of the steps that I might take:
The $4\times 4$ matrix $H$ can always be written in the form
$$
H=a\mathbb{I}\otimes\mathbb{I}+\underline{n}_1\cdot\underline{\sigma}\otimes\mathbb{I}+\mathbb{I}\otimes\underline{n}_2\cdot\underline{\sigma}+\underline{\sigma}\cdot M\cdot\underline{\sigma}
$$
where $M$ is a real $3\times 3$ matrix. Remember that if you want to find the coefficient of a particular term, then you can calculate
$$
\text{Tr}(H(\sigma_i\otimes\sigma_j))/4
$$
for $i,j\in\{0,1,2,3\}$.
Now, if you implement the unitary $U_i$ on qubit $i$, Pauli matrices change as
$$
\sigma_j\mapsto \underline{R}_j^{(i)}\cdot\underline{\sigma},
$$
and, if you work it through, you find that $M$ updates to $R^{(1)}\cdot M\cdot {R^{(2)}}^T$. So, you can choose the $R$ matrices to be the matrices that yield the singular values of $M$. In that way, you only ever have to deal with a matrix of the form
$$
H= a\mathbb{I}\otimes\mathbb{I}+\underline{m}_1\cdot\underline{\sigma}\otimes\mathbb{I}+\mathbb{I}\otimes\underline{m}_2\cdot\underline{\sigma}+n_1X\otimes X+n_2Y\otimes Y+n_3Z\otimes Z.
$$
If you get lucky and the local fields are in the Z-direction only, this matrix divides into two $2\times 2$ matrices spanned by $\{|00\rangle,|11\rangle\}$ and $\{|01\rangle,|10\rangle\}$ respectively.
The other option is that if you get really lucky (e.g. in this question!), the non-zero terms in the decomposition all mutually commute. Then you can analytically diagonalise $H$ in a much easier manner.
However, all of that is far more work that just throwing the $4\times 4$ matrix into the computer and asking for the eigenvalues. After all, I've already had to require diagonalization of a $3\times 3$ matrix. (Where it is more useful is if you have a translation invariant Hamiltonian of many qubits.)