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.
Best Answer
Diagonalizing the Hamiltonian means you want to bring it into the form $H=\omega b^\dagger b$, and it is pretty obvious that $b$ should be a linear combination of $a$ and $a^\dagger$, and $b$ should satisfy the canonical commutation of annihilation operators, namely $[b,b^\dagger]=1, [b,b]=0$.
Now let's write $b=ua+va^\dagger$ (this is called the Bogoliubov transformation, by the way). The condition $[b,b^\dagger]=1$ leads to $|u|^2-|v|^2=1$. Let us expand out $b^\dagger b$:
$$ b^\dagger b= |u|^2 a^\dagger a+ |v|^2 a a^\dagger + u^*v a^\dagger a^\dagger + uv^* aa. $$
Therefore
$$ \omega(|u|^2+|v|^2)=E+\Delta, \omega u^*v = \frac{1}{2}\Delta. $$
Together with $|u|^2-|v|^2=1$, we have three equations for three variables ($u, v, \omega$). In fact, in this case one can safely assume $u$ and $v$ are both real. The rest is just algebra.