Is there a way to express the general $4 \times 4$ Hamiltonian in some block diagonal form of $2 \times 2$ matrices that I can solve, knowing the exact solution of $2\times 2$?
This is necessary for the treatment I am going to perform later, as I cannot go for general quartics solutions of this system. I need to see the action of the Hamiltonian to subspaces and the interaction between them.
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.)
If you want to diagonalize a general 4x4 Hamiltonian, you cannot just diagonalize the four 2x2 blocks and piece together what you get. You need to diagonalize the entire 4x4 matrix all at once.
If the 4x4 matrix is already block diagonal, then of course you can diagonalize each 2x2 block separately, but not in the general case where all 16 elements of the 4x4 matrix can be arbitrary.
