Take the Hamiltonian you write at end of your question as the example(where your statement about the eigenvector is not right), the procedure to obtain the eigenvalue and eigenvectors of other Hamiltonian is nearly the same.
$$H=\sum_{ij}A_{ij}a_i^\dagger a_j$$
Writing in matrix form you have(suppose it is a $n\times n$ matrix):
$$H=\alpha^\dagger A \alpha$$
with:
$$ \alpha^\dagger=(a_1^\dagger, a_2^\dagger,\dotsb,a_n^\dagger )$$
Now you are ready to diagonalize Matrix analytically or numerically:
$$A=X^\dagger D X$$
Where $D=[E_1,E_2,\dotsb,E_n]$ is the eigenvalue you want, and $X$ is just a Unitary matrix that diagonalize $A$. Indeed, it is the normalized mathematical eigenvectors of $A$ write in column way(which is same in Fortran if you do it numerically).
Now substitute this to the original Hamiltonian:
$$H=\alpha^\dagger A \alpha=\alpha^\dagger X^\dagger DX\alpha=\beta^\dagger D \beta$$
where $\beta$ is $\beta =X\alpha=(\beta_1,\beta_2,\dotsb,\beta_n)^T$. Because this is a unitary transformation, you can easily check that the commute or anti-commute relation still holds for $\beta_i,\beta_2,\dotsb$.
Now you can consider the system is consisting of non-interacting qusi-particles $\beta_i,\beta_2,\dotsb$, then the eigenstate of the system is easily obtained:
$$|\psi \rangle=(\beta_1^\dagger)^{n_1}(\beta_2^\dagger)^{n_2}\dotsb(\beta_n^\dagger)^{n_n}|0\rangle$$
$n_i=0,1$ for fermions and $n_i=0,1,2,3\dotsb$ for bosons.
Also, see this thread with the very similar question.
You can try just working out the matrix elements between the different localised states, i.e. $a_1^{\dagger}\ldots a_m^{\dagger}|0\rangle$ etc. This may be rather tedious and require you to carefully keep track of phases due to the fermionic anti-symmetry. Probably a better way is to use the Jordan-Wigner transformation, which gives you an explicit matrix representation of the fermion operators. In 2D this representation of the Hamiltonian looks ugly, but it is still perfectly useful. You will be restricted to a fairly small lattice size, i.e. 10-20 sites depending on your computer and the efficiency of your code. It will be more efficient if you use the fermion number conservation to diagonalise within each number sector separately.
In the particularly simple case of a single particle on the lattice, the problem becomes rather easy, since exchange statistics are no longer relevant. In this case, you make the replacement
$$a_{i,j}^{\dagger}a_{m,n} \to |i,j\rangle\langle m,n|,$$
where
$$|i,j\rangle = a_{i,j}^{\dagger}|0\rangle.$$
Best Answer
The energies and states you get are correct. I will point out something else that may be useful. I will use the following basis: $$|00> = |0>;\quad|01> = d^\dagger|0>; |10> = f^\dagger|0>;\quad |11>=f^\dagger d^\dagger |0>.$$ The Hamiltonian at hand conserves parity, so you can write the Fock-space matrix in a more convenient way: $$H=\begin{array}{c|cc|cc} \,&|00>&|11>& |01> & |10> \\ \hline <00| & -(\epsilon_m + \epsilon_d/2) & t & 0 & 0\\ <11| & t & (\epsilon_m + \epsilon_d/2) & 0 & 0 \\ \hline <01|&0 & 0& -(\epsilon_m-\epsilon_d/2) & t \\ <10| & 0 & 0 & t & \epsilon_m-\epsilon_d/2 \end{array},$$ which is block-diagonal, so you only have to deal with subspaces of size $2\times 2$. Let's label the parities of the subspaces as $P=0$ and $P=1$, respectively, as $P=(n_f+n_d)\text{mod}2$. The eigenvalues are $$E_{0,\pm}=\pm\frac{\sqrt{(\epsilon_d+2\epsilon_m)^2+4t^2}}{2};\quad E_{1,\pm}=\pm\frac{\sqrt{(\epsilon_d-2\epsilon_m)^2+4t^2}}{2},$$ where the index $0$ or $1$ refers to the parity. You already got the normalized eigenstates, and all you need to do is reorganize them to put them in this particular ordering. To make things quicker, I'll simply write them as follows: $$|0\pm>=\alpha_0^{\pm}|00>+\beta_0^{\pm}|11>;\quad |1\pm>=\alpha_1^{\pm}|01>+\beta_1^{\pm}|10>,$$ with $\alpha_P^{\pm}$ and $\beta_{P}^{\pm}$ the coefficients you found. You can associate to these the operators $A_{0,\pm}^\dagger|0> = |0\pm>$ and $A_{1,\pm}^\dagger|0>=|1\pm>$.
Now, the operators $A_{1,\pm}$ have fermion commutation relations, because they are single-particle operators: $$\{A_{1,\sigma}^\dagger,A_{1,\sigma'}\}=\{\alpha_{1,\sigma}d^\dagger+\beta_{1,\sigma}f^\dagger,\alpha_{1,\sigma'}^*d+\beta_{1,\sigma'}^*f\}=\alpha_{1,\sigma}\alpha_{1,\sigma'}^*+\beta_{1,\sigma}\beta_{1,\sigma'}^*=\delta_{\sigma,\sigma'},$$ by definition (the coefficients meet this condition, because they are the coefficients of mutually orthonormal vectors). For the $P=0$ subspace, what you have are actually bosons, but not normal bosons, because they mix different parts of Fock space ($0$ particle and $2$ particle). As a consequence, you cannot get normal boson commutation relations: $$\begin{split}[A_{0,\sigma}^\dagger,A_{0,\sigma'}]=&[\alpha_{0,\sigma}+\beta_{0,\sigma}f^\dagger d^\dagger,\alpha_{0,\sigma'}^*+\beta_{0,\sigma'}^*d f]\\ =& \beta_{0\sigma}\beta_{0\sigma'}^*\left(d^\dagger d + f^\dagger f -1 \right).\end{split}$$ What you can do, on the other hand, is define operators of the form $B=\gamma d + \delta f^\dagger$ and $C = \lambda d^\dagger + \mu f$, such that you can make products of the form $$B A_{1,\sigma}^\dagger = \gamma \alpha_{1,\sigma} (1-n_d) + \gamma \beta_{1,\sigma} d f^\dagger + \delta \alpha_{1,\sigma} f^\dagger d^\dagger + \delta \beta_{1,\sigma} (1-n_f).$$ Applied to the vacuum, the terms $d f^\dagger$, $n_f$ and $n_d$ vanish, and you get $$B A^\dagger_{1,\sigma}|0> = (\gamma \alpha_{1\sigma} + \delta \beta_{1\sigma})|00> + \delta \alpha_{1\sigma}|11>.$$ The coefficients $\gamma$ and $\delta$ can be fixed such that $$(\gamma \alpha_{1\sigma} + \delta \beta_{1\sigma})=\alpha_{0,\sigma'}\quad;\quad \delta \alpha_{1\sigma} = \beta_{0,\sigma'}.$$ The weird fermion $B$ satisfies $\{B,B^\dagger\}=|\gamma|^2+|\delta|^2$.
I'll leave the rest for you to work out. BTW, the Hamiltonian you are solving represents the coupling between a Majorana fermion $\gamma_1$, and an ordinary fermion $(d^\dagger, d)$, in the form $H=\epsilon_d d^\dagger d + \gamma_1 (d^\dagger - d)$. The Majorana fermion can be written in terms of the auxiliary fermion $(f^\dagger,f)$ as $\gamma_1 = f^\dagger + f$. You can define an independent Majorana $i\gamma_2 = f^\dagger - f$, and together they satisfy the Clifford algebra $\{\gamma_i, \gamma_j \}=2\delta_{ij}$. If you have further questions, or want to learn more, you can look up papers about Majoranas coupled to normal fermions. I hope this was of some use to you.