The working matrix $(n+1)\times (n+1)$:
$$
\mathbf{H} =
\begin{bmatrix}
\epsilon_0 & U_1 & U_2 & U_3 & ... &U_{n-1}&U_n\\
U_1 & \epsilon_1 &0&0&...&0&0\\
U_2& 0 & \epsilon_2 &0&...&0&0\\
U_3& 0 & 0 & \epsilon_3 &...&0&0\\
...&...&...&...&...&....&...\\
U_{n-1} &0&0&0&...&\epsilon_{n-1}&0\\
U_n &0&0&0&...&0&\epsilon_n
\end{bmatrix}
$$
Assum that all $U_i \ne = 0$, and are real numbers. The eigen equation for this matrix:
$$ \tag{1}
det\left( \mathbf{H} - \lambda \mathbf{I} \right) = 0
$$
Since $\mathbf{H}$ is a symmetric real matrix, it guarantees that there exists $n+1$ real eigen values and the corresponding $n+1$ orthogonal eigen vectors.
A closed form of the eigen equation can be obtained by Gauss-Jordan elimination. First, add to the first row, the $n$th row mutiplied by $-\frac{U_n}{\epsilon_n -\lambda}$. This will leave $H_{0n}=0$,
and adds a term to $H_{00}$, $H_{00} = (\epsilon_0 - \lambda) -\frac{U_n^2}{\epsilon_n -\lambda}$.
Repeat this process for $i = n-1, n-2, ... 1$. After $n$ eliminations, the upper triangular elements of $\mathbf{H} - \lambda \mathbf{I}$ becomes all zeros, and the element $H_{00}$ collects $n$ more terms from the elimination process:
$$ \tag{2}
H_{00} = (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda}.
$$
and the determinant is the product of all diagonal elements:
$$
det\left(\mathbf{H} - \lambda \mathbf{I} \right) = \left\{ (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda} \right\} \prod_{i=1}^n (\epsilon_i - \lambda).
$$
$$ \tag{3}
= \prod_{i=0}^n (\epsilon_i - \lambda) - \sum_{i=1}^n U^2_i \prod_{j=1, j\ne i}^n (\epsilon_j - \lambda).
$$
Case 1. Degenerate diaginal elements
In case that there are some energies are equal, saying 3-fold degeneracy: $\epsilon_a = \epsilon_b = \epsilon_c$
for $a$, $b$, $c$ are arbitrary three integers in $\{1,2,3, ..., n\}$.
Under this case, a term $(\epsilon_a - \lambda )^2$ can be factored out of the Eq. (3).
The energy $\epsilon_a$ is therefore an eigen value of degeneracy $3-1=2$. The corresponding
eigen vectors can be calculated from the $4\times 4$ submatrix:
$$ \tag{4}
\begin{bmatrix}
\epsilon_0-\epsilon_a & U_a & U_b & U_c\\
U_a & 0 & 0 & 0\\
U_b & 0 & 0 & 0\\
U_c & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
x_0\\
x_a\\
x_b\\
x_c
\end{bmatrix}
=
\begin{bmatrix}
0\\
0\\
0\\
0
\end{bmatrix}
$$
The equation for eigen vectors:
$$ x_0 = 0;$$
$$ x_a U_a + x_b U_b + x_c U_c = 0. $$
Two orthonomal eigen vectors can be constructed from these two conditions. It would be very interesting to observe how these two domension block-diaginalized and decoupled from the rest $(n-1)\times (n-1)$ dimensions, and how the coupling elements is changed for the remaining state of the 3 degenerate levels.
We conclude that for $p$-fold degeneracy energy $\epsilon_a$, ther existed $(p-1)$ degenerate eigen vectors with eigen value $\lambda = \epsilon_a$. This $(p-1)$ basis vectors span an $(p-1)$ subspace. The motion of the subspace is decoupled from the rest of the system. One example realized this matrix is a coupled oscillator system, which all the oscillators are coupled only to the $0$th oscillator (see Figure). The dynamical matrix for the normal mode will resemble this matrix form. Then, we can select a subgroup of oscillation modes, that they are coupled only within the subgroup.
Another worthy mentioned situation is that all $\epsilon_i = \epsilon_1$ for all $i=1, 2, 3,..,n$. Then we immediately has $n-1$-fold eigen values $\lambda = \epsilon_1$. Their eigen vectors span the $(n-1)\times (n-1)$ subspace.
$$ \vec{v}_\lambda = [x_0, x_1, x_2, ..., x_{n-1}, x_n] $$
Conditions for eigen vector:
$$ x_0 = 0;$$
$$ \sum_{i=1}^n x_i U_i = 0. $$
Since we have $(n+1)\times (n+1)$ symmetric matrix, there are two more eigen values. From the determinant, after factoring out $(\epsilon_1-\lambda)^{n-1}$, we find that the last two eigen values are the root of the equation:
$$
(\epsilon_0-\lambda) (\epsilon_1-\lambda) - \sum_{i=1}^n U_n^2 = 0.
$$
Which is two-level matrix with coupling strength $U =\sqrt{ \sum_{i=1}^n U_n^2} $. The eigen values are readily given:
$$
\lambda_\pm = \frac{\epsilon_0+\epsilon_1}{2} \pm \sqrt{ \frac{(\epsilon_0-\epsilon_1)^2}{2} + U^2}
$$
The total $n+1$ eigen values consist of $(n-1)$-fold of $\epsilon_1$ in the middle and two $ \lambda_\pm$ in higher and lower values.
Non-degenrate case
we may now assume that all $\epsilon_i$ are different. All eigen values will then solely determined by zeros of the element $H_{00}$ in Eq.(2)':
{2}
$$
H_{00} = (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda} = 0.
$$
$$
= (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{(\epsilon_i-\epsilon_0 ) -( \lambda -\epsilon_0 ) } = 0.
$$
Setting $\epsilon_0 = 0$ without lossing generality:
$$
= -\lambda - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda } = 0.
$$
For better understanding of this equation, lets assum $\epsilon_i = i $. We can now examine the value of $H_{00}(\lambda)$ by scanning $\lambda$ from $-\infty$ to $\infty$.
The value of $H_{00}(\lambda)$ is positive when $\lambda \ll 0 $, but the term of summation is negative. Just before $\lambda \to 0^-$, the smaller negative summation catches up with $-\lambda$. It is the smallest eigen value (the $\lambda_-$ in the total degenerate case.)
Then the $H_{00}$ approach $-\infty$ at $\epsilon_1^-$ due to divergence in
$-\frac{U_1^2}{\epsilon_1 - \lambda}$. After $\lambda = \epsilon_1^+$, the function emerges with large positive value. And approach $-\infty$ again at $\epsilon_2^-$, behaves like the $\cot$ function. Therefore, it locates first eigen value at $\lambda_0 <0$, and $n-1$ eigen values at $\epsilon_i < \lambda_i < \epsilon_{i+1}$, a final eigen value locates
at $\lambda_n > \epsilon_n $, because when
$\lambda \to \infty$, $H_{00} \to -\infty$, as shown in Fig.(2).
We can see this figure topologically, as $\epsilon_1$,..$\epsilon_n$ are drawn closer, and finally into one energy $\epsilon_1$. All the levels in between, are forced collapse into a single degenerate level, but the first and the last, stay away from this squeezing process.
pertubation for small $U_i$'s
In case that $U_i \ll |\epsilon_a - \epsilon_b|$ for any two diagonal element, $\epsilon_a$, and
$\epsilon_b$. The eigen value and eigen vector can be obtained by a approximation method similar
to the second order perturbation in quantum mechanics:
$$
\lambda_0 = \sum_{i=0}^n \frac{U_i^2}{\epsilon_i}.
$$
Where is replace the $\lambda$ in the summation by $\epsilon_0 = 0$.
For $i \ge 1$
$$
\lambda_i = \epsilon_i + \frac{U_i^2}{\epsilon_i - \epsilon_0} =
\epsilon_i + \frac{U_i^2}{\epsilon_i}.
$$
Where we replace the $\lambda$ in the summation by $\epsilon_i$, and made expansion for small $U_i$s.
Eigen vectors and tranformation
Each eigen value renders an eigen vector $\vec{v}_\lambda = \left[ x_0, x_1, x_2, ...,x_n \right] $:
$$
\vec{v}_\lambda = A_1
\begin{bmatrix}
1, \frac{U_1}{\lambda - \epsilon_1}, \frac{U_2}{\lambda - \epsilon_2}, ..., \frac{U_n}{\lambda - \epsilon_n}
\end{bmatrix}
$$
Leave $A_1$ as the normalization constant. There are $n$ ortho-normal vectors, providing the bases for the orthogonal transformation $\mathbf{R}$ to the diagonalisation of the Hamiltonian.
$$
\mathbf{R} =
\begin{bmatrix}
A_0 & A_1 & A_2 & A_3& ... & A_n\\
\frac{A_0 U_1}{\lambda_0 - \epsilon_1} &\frac{A_1 U_1}{\lambda_1 - \epsilon_1} & \frac{A_2 U_2}{\lambda_2 - \epsilon_2} &\frac{A_3 U_1}{\lambda_3 - \epsilon_1} & ... & \frac{A_{n} U_1}{\lambda_{n} - \epsilon_1}\\
\frac{A_0 U_2}{\lambda_0 - \epsilon_1} &\frac{A_1 U_2}{\lambda_1 - \epsilon_2} & \frac{A_2 U_2}{\lambda_2 - \epsilon_1} &\frac{A_3 U_2}{\lambda_3 - \epsilon_2} & ... & \frac{A_{n} U_2}{\lambda_{n} - \epsilon_2}\\
....&... & ... & ... & ... \\
\frac{A_0 U_{n}}{\lambda_0 - \epsilon_1} & \frac{A_1 U_{n}}{\lambda_1 - \epsilon_{n}} & \frac{A_2 U_{n}}{\lambda_2 - \epsilon_{n}} & \frac{A_3 U_{n}}{\lambda_3 - \epsilon_{n}} &... & \frac{A_{n} U_{n}}{\lambda_{n} - \epsilon_{n}}
\end{bmatrix}.
$$
Since the $\mathbf{H}$ is Hermitian, the $n+1$ eigen vectors are mutually orthogonal:
$$
\mathbf{R} \mathbf{R}^T =\mathbf{R}^T \mathbf{R} = \mathbf{I}
$$
$$
\mathbf{R}^T \mathbf{H}\mathbf{R} =\mathbf{D}
$$
Where $\mathbf{D}$ is the diagonal matrix with all $\lambda_i$ as the elements.
Lets consider a matrix equation:
$$
\mathbf{H} \vec{v} = \vec{w}
$$
Transform it to the eigen coordinate:
$$
\mathbf{R}^T \mathbf{H} \left(\mathbf{R} \mathbf{R}^T \right) \vec{v} = \mathbf{R}^T\vec{w}
$$
$$
\mathbf{D} \left( \mathbf{R}^T \vec{v} \right) = \left(\mathbf{R}^T\vec{w}\right)
$$
$$
\mathbf{D} \vec{v}' = \vec{w}'
$$
In $\vec{v}'$ and $ \vec{w}'$ primed coordinate, the matrix $\mathbf{D}$ is diagonalized, thus each components of $\vec{v}'$ and $ \vec{w}'$ is decoupled.
The transformation
$$
\mathbf{R}^T \vec{w} = \vec{w}' ;
$$
and
$$
\mathbf{R} \vec{w}' = \vec{w} ;
$$
