Your Hilbert space is finite-dimensional (specifically, 4-dimensional). It means that you need no fancy way of deriving the spectrum, just proceed with the standard procedure. First, write your Hamiltonian in the matrix form:
$$ H = \left(
\begin{array}{cccc}
\epsilon & -\Delta & 0 & 0 \\
-\Delta & 0 & 0 & 0 \\
0 & 0 & \epsilon & -\Delta \\
0 & 0 & -\Delta & 0 \\
\end{array}
\right). $$
You can immediately see that your system is actually a direct sum of two independent systems with hamiltonians
$$ h = \left(
\begin{array}{cc}
\epsilon & -\Delta \\
-\Delta & 0 \\
\end{array}
\right). $$
But let me proceed as if I don't know that.
You can solve the characteristic equation $\det \left( H - \lambda \cdot 1_{4 \times 4} \right) = 0$ with respect to $\lambda$ which would give you the eigenvalues (energy spectrum) of $H$:
$$ \left\{\frac{1}{2} \left(\epsilon -\sqrt{4 \Delta ^2+\epsilon
^2}\right),\frac{1}{2} \left(\sqrt{4 \Delta ^2+\epsilon ^2}+\epsilon
\right)\right\}. $$
Both eigenstates are degenerate of order 2 (because they are roots of order 2 of the characteristic equation).
The corresponding eigenstates are:
$$ \left\{\left(
\begin{array}{c}
0 \\
0 \\
\frac{\sqrt{4 \Delta ^2+\epsilon ^2}-\epsilon }{2 \Delta } \\
1 \\
\end{array}
\right),\left(
\begin{array}{c}
\frac{\sqrt{4 \Delta ^2+\epsilon ^2}-\epsilon }{2 \Delta } \\
1 \\
0 \\
0 \\
\end{array}
\right)\right\} $$
(for the first eigenvalue) and
$$ \left\{\left(
\begin{array}{c}
0 \\
0 \\
-\frac{\epsilon +\sqrt{4 \Delta ^2+\epsilon ^2}}{2 \Delta } \\
1 \\
\end{array}
\right),\left(
\begin{array}{c}
-\frac{\epsilon +\sqrt{4 \Delta ^2+\epsilon ^2}}{2 \Delta } \\
1 \\
0 \\
0 \\
\end{array}
\right)\right\} $$
for the second one.
No perturbation expansion is needed! You can build an arbitrary state time evolution out of the energy eigenstates via the general formula:
$$ \Psi = \sum e^{-i E_a t / \hbar} \psi_a $$
Note that $\epsilon$ and $\Delta$ have to be real in order for $H$ to be hermitian.