The example you find has a mistake. When you diagonalize the matrix $\mathcal{H}$, you are applying some unitary transformation $U$, which is a $4 \times 4$ matrix, to the vector $(c_1,c^{\dagger}_2,c_2,c^{\dagger}_1)^T$. Let the transformed vector be $(\tilde{c}_1,\tilde{c}^{\dagger}_2,\tilde{c}_2,\tilde{c}^{\dagger}_1)^T$. We must have $\tilde{c}^{\dagger}_1=(\tilde{c}_1)^{\dagger}$ and $\tilde{c}^{\dagger}_2=(\tilde{c}_2)^{\dagger}$, but a $4 \times 4$ matrix $U$ gives extra degrees of freedoms to the transformation, and the relation generally no longer holds. Therefore, you will have different eigenvalues for the same Hamiltonian as invalid results.

In your previous edition of the post, you mentioned a link about [Jordan-Wigner transformation][1] and its use in $1$-dimensional quantum Ising model. The Hamiltonian after Jordan-Wigner transformation is
$$H=\sum_{k}{\big(-c^{\dagger}_kc_k(2\cos{k})-(c^{\dagger}_kc^{\dagger}_{-k}e^{ik}+c_{-k}c_ke^{ik})+2h_zc^{\dagger}_kc_k\big)}$$
You argue that we can write the Hamiltonian as
$$H=\sum_{k}
\begin{pmatrix} c^{\dagger}_k & c_{-k} \end{pmatrix}
\begin{pmatrix}
-2\cos{k}+2h_z & -e^{ik} \\
-e^{-ik} & 0
\end{pmatrix}
\begin{pmatrix} c_k \\ c^{\dagger}_{-k} \end{pmatrix}$$
This is true, but diagonalizing the $2x2$ matrix in the above equation does not give you superposition of $c_k$ and $c^{\dagger}_{-k}$ which can diagonalize the Hamiltonian. There is the other $-k$ term as contribution in the sum, meaning
\begin{align}
H & = \sum_{k>0}
\begin{pmatrix} c^{\dagger}_k & c_{-k} \end{pmatrix}
\begin{pmatrix}
-2\cos{k}+2h_z & -e^{ik} \\
-e^{-ik} & 0
\end{pmatrix}
\begin{pmatrix} c_k \\ c^{\dagger}_{-k} \end{pmatrix} \\
& + \sum_{-k,k \geq 0}
\begin{pmatrix} c^{\dagger}_{-k} & c_k \end{pmatrix}
\begin{pmatrix}
-2\cos{k}+2h_z & -e^{-ik} \\
-e^{ik} & 0
\end{pmatrix}
\begin{pmatrix} c_{-k} \\ c^{\dagger}_k \end{pmatrix}
\end{align}
Therefore, you can see the diagonalization of $\begin{pmatrix}
-2\cos{k}+2h_z & -e^{ik} \\
-e^{-ik} & 0
\end{pmatrix}$ only diagonalizes the first part of the Hamiltonian. As written in the link, the correct way to do this is to sum up the two parts, and have
$$H=\sum_{k}
\begin{pmatrix} c^{\dagger}_k & c_{-k} \end{pmatrix}
\begin{pmatrix}
-\cos{k}+h_z & -i\sin{k} \\
i\sin{k} & \cos{k}-h_z
\end{pmatrix}
\begin{pmatrix} c_k \\ c^{\dagger}_{-k} \end{pmatrix}$$
The diagonalization of the matrix gives you the superposition of $c_k$ and $c^{\dagger}$ which diagonalizes the Hamiltonian.


 


  [1]: https://theory.leeds.ac.uk/interaction-distance/applications/ising/map-to-free/