Assume the tridiagonal matrix $T$ is in this form:
$$
T = \begin{bmatrix}
a & c & & & &\\
b & a & c & & &\\
& b & a & c & &\\
& & &\ddots & &\\
& & & b & a & c\\
& & & & b & a\\
\end{bmatrix}
$$
we must show that its eigenvalues are of the form
$$a + 2 \sqrt{bc} \, \cos \left( \frac{k \pi}{n+1} \right)$$
where $$a=qh^2−1, ~~ b=1- \frac{ph}{2}, ~~ c =1+\frac{ph}{2} , ~~q \leq 0.$$
Best Answer
I am not sure that this s entirely rigorous (!!!) but here's an idea.
Let $T_n$ be your tridiagonal matrix of order $n$, and let $S_n=T_n-\mathbb{I}\sigma$. Let $d_n$ be the determinant of $S_n$. Solving $d_n=0$ gives the desired eigenvalues $\sigma_1,\dots,\sigma_n$. Developing $d_n$ with Laplace's rule and letting $a'=a-\sigma$, you have the recurrence relation $d_{n+1}=a'\cdot d_{n}-bc\cdot d_{n-1}$. You can assume $d_0=1$ and $d_1=a'$. You can now set up the following system: \begin{equation} \left(\begin{array}{c} d_{n+1}\\ d_{n} \end{array}\right) = \left( \begin{array}{cc} a' & -bc \\ 1 & 0 \end{array} \right) \left(\begin{array}{c} d_n\\ d_{n-1} \end{array}\right) \end{equation} Let $v_{n}=(d_{n+1},d_n)$ and $A$ the above matrix. So you have \begin{equation} v_n=Av_{n-1}=A^2v_{n-2}=\cdots=A^nv_0 \end{equation} where $v_0=(d_1,d_0)=(a',1)$. Your problem now is getting $A^n$. You can use the Cayley-Hamilton theorem to write \begin{equation} A^n=\alpha A+\beta\mathbb{I} \end{equation} To find $\alpha$ and $\beta$ you need the eigenvalues of $A$, let them be $\lambda_1$ and $\lambda_2$. Then, since the above equation is obtained from the characteristic polynomial, it must be satisfied also by the eigenvalues, so the following system holds: \begin{equation} \lambda_1^n=\alpha \lambda_1+\beta\\ \lambda_2^n=\alpha \lambda_2+\beta \end{equation} If the two eigenvalues are different, you get $\alpha=\frac{\lambda_1^n-\lambda_2^n}{\lambda_1-\lambda_2}$ and $\beta=-\lambda_1\lambda_2\frac{\lambda_1^{n-1}-\lambda_2^{n-1}}{\lambda_1-\lambda_2}$, so \begin{equation} v_n=\alpha Av_0+\beta v_0\quad\Longrightarrow\quad d_n=\alpha a'+\beta = \frac{\lambda_1^n-\lambda_2^n}{\lambda_1-\lambda_2}a'-\lambda_1\lambda_2\frac{\lambda_1^{n-1}-\lambda_2^{n-1}}{\lambda_1-\lambda_2} \end{equation} Simplifying, \begin{equation} d_n=\frac{a'\lambda_1^n-a'\lambda_2^n-\lambda_2\lambda_1^n+\lambda_1\lambda_2^n}{\lambda_1-\lambda_2}=\frac{\lambda_1^n(a'-\lambda_2)-\lambda_2^n(a'-\lambda_1)}{\lambda_1-\lambda_2} \end{equation} Computing the eigenvalues you get $\lambda_1 = \frac{1}{2}(a'+\sqrt{a'^2-4bc})$ and $\lambda_2=\frac{1}{2}(a'-\sqrt{a'^2-4bc})$ so we have $\lambda_1+\lambda_2=a'$, and we can write \begin{equation} d_n=\frac{\lambda_1^{n+1}-\lambda_2^{n+1}}{\lambda_1-\lambda_2} \end{equation} Since the eigenvalues were different, we can suppose without loss of generality that $\lambda_2\neq 0$, so: \begin{equation} d_n=\lambda_2^n\frac{(\lambda_1/\lambda_2)^{n+1}-1}{(\lambda_1/\lambda_2)-1} \end{equation} Now - the original problem was finding the eigenvalues of the big matrix. That means that we must search for the solutions of $d_n=0$, i.e. \begin{equation} (\lambda_1/\lambda_2)^{n+1}=1 \end{equation} Luckily we have invented complex numbers, so \begin{equation} \lambda_1/\lambda_2=e^{i\frac{2k\pi}{n+1}}\qquad\Longrightarrow\qquad\lambda_1=\lambda_2e^{i\frac{2k\pi}{n+1}}\qquad \text{for }k\in[1\cdots n] \end{equation} Note that we cannot let $k=0$ because we supposed $\lambda_1\neq\lambda_2$. Calling $z=e^{i\frac{2k\pi}{n+1}}$ and substituting the expression for $\lambda_1$ and $\lambda_2$ we get \begin{equation} a'+\sqrt{a'^2-4bc}=(a'-\sqrt{a'^2-4bc})z \end{equation} or \begin{equation} a'\cdot(z-1)=\sqrt{a'^2-4bc}\cdot(z+1) \end{equation} Squaring, \begin{equation} a'^2(z-1)^2=a'^2(z+1)^2-4bc(z+1)^2\Longrightarrow a'^2=\frac{(z+1)^2}{z}bc \end{equation} So \begin{equation} a'=\frac{z+1}{\sqrt{z}}\sqrt{bc}=\frac{e^{i\frac{2k\pi}{n+1}}+1}{e^{i\frac{k\pi}{n+1}}}\sqrt{bc} \end{equation} Multiplying by $1=e^{-i\frac{k\pi}{n+1}}/e^{-i\frac{k\pi}{n+1}}$ we get \begin{equation} a'=(e^{i\frac{k\pi}{n+1}}+e^{-i\frac{k\pi}{n+1}})\sqrt{bc} \end{equation} The complex exponential is actually real since $e^{i\theta}+e^{-i\theta}=2\cos\theta$. We had also $a'=a-\sigma$, so the result follows: \begin{equation} \sigma=a-2\cos\left(\frac{k\pi}{n+1}\right)\sqrt{bc}=a+2\cos\left(\pi-\frac{k\pi}{n+1}\right)\sqrt{bc}=a+2\cos\left(\frac{h\pi}{n+1}\right)\sqrt{bc} \end{equation} with $h\in [1,n]$.
EDIT: what happens if $\lambda_1=\lambda_2=\lambda$? In this case, we have a single eigenvalue so the characteristic polynomial of $A$ has the form $p(x)=(x-\lambda)^2$. We see that $p'(x)=2(x-\lambda)x$ so both $p(\lambda)$ and $p'(\lambda)$ are zero. Then, we can consider the system \begin{equation} \lambda^n=\alpha \lambda+\beta\\ n\lambda^{n-1}=\alpha \end{equation} so $\alpha=n\lambda^{n-1}$ and $\beta=(1-n)\lambda^n$, with $\lambda=a'/2$.