First of all you should recall that Schroedinger equation is an Eigenvalue equation. If you are unfamiliar with eigenvalue equations you should consult any math book or course as soon as possible.
Answer 1 (my apologies, I will use my own notation, as this is mainly copy-paste from my old notes):
First define constants
\begin{equation}
x_0 = \sqrt{\frac{\hbar}{m\omega}} ,
\end{equation}
\begin{equation}
p_0 = \frac{\hbar}{x_0} = \sqrt{\hbar m \omega} ,
\end{equation}
and dimensionless operators
\begin{equation}
\hat{X} = \frac{1}{x_0} \hat{x} ,
\end{equation}
and
\begin{equation}
\hat{P} = \frac{1}{p_0} \hat{p} .
\end{equation}
Their commutation relation then is
\begin{equation}
\left[ \hat{X} , \hat{P} \right] = \left[ \frac{1}{x_0}\hat{x} , \frac{1}{p_0}\hat{p} \right] = \frac{1}{x_0 p_0} \left(\hat{x} \hat{p} - \hat{p} \hat{x} \right) = \frac{1}{x_0 p_0} \left[ \hat{x} , \hat{p} \right] = \frac{i\hbar}{x_0 p_0} = i ,
\end{equation}
as
\begin{equation}
x_0 p_0 = \sqrt{\frac{\hbar}{m\omega}} \sqrt{\hbar m\omega} = \hbar .
\end{equation}
Now write Hamiltonian in terms of $\hat{X}$ and $\hat{P}$. Start with
\begin{equation}
\hat{H} = \frac{p_0 ^2}{2m} \hat{P} ^2 + \frac{1}{2} m\omega^2 x_0 ^2 \hat{X}^2 .
\end{equation}
Notice that
\begin{equation}
p_0 ^2 = \hbar m \omega
\end{equation}
and
\begin{equation}
x_0 ^2 = \frac{\hbar}{m \omega} ,
\end{equation}
hence
\begin{equation}
\hat{H} = \frac{\hbar \omega}{2} \hat{P}^2 + \frac{\hbar \omega}{2} \hat{X}^2 = \frac{\hbar \omega}{2} \left(\hat{X}^2 + \hat{P}^2 \right).
\end{equation}
Up to the commutation relation we can write
\begin{equation}
\left(X^2 + P^2 \right) = \left(X - iP \right) \left(X + iP \right) .
\end{equation}
On the other hand, for operators this is not quite allowed, as
\begin{align}
\left(\hat{X} - i\hat{P} \right) \left(\hat{X} + i\hat{P} \right) &= \hat{X}^2 + i\hat{X}\hat{P} - i\hat{P}\hat{X} + \hat{P}^2 \nonumber \\ &= \hat{X}^2 +i \left(\hat{X}\hat{P} - \hat{P}\hat{X}\right) + \hat{P}^2 \\ &= \hat{X}^2 + i \left[ \hat{X} , \hat{P} \right] + \hat{P}^2 = \hat{X}^2 + \hat{P}^2 - 1 \nonumber ,
\end{align}
so one has
\begin{equation}
\left(\hat{X}^2 + \hat{P}^2 \right) = \left(\hat{X} - i\hat{P} \right) \left(\hat{X} + i\hat{P} \right) + 1 .
\end{equation}
Now we can define
\begin{equation}
\hat{a} = \frac{1}{\sqrt{2}} \left(\hat{X} + i\hat{P} \right) ,
\end{equation}
and
\begin{equation}
\hat{a}^\dagger = \frac{1}{\sqrt{2}} \left(\hat{X} - i\hat{P} \right),
\end{equation}
calling this creation operator and $\hat{a}$ - annihilation operator. Notice that we can now express Hamiltonian in terms of creation and annihilation operators:
\begin{equation}
\hat{H} = \frac{\hbar\omega}{2} \left(\sqrt{2} \hat{a} ^\dagger \sqrt{2} \hat{a} + 1 \right) = \hbar \omega \left(\hat{a} ^\dagger \hat{a} + \frac{1}{2} \right) .
\end{equation}
But we can also define the number operator, $\hat{N} = \hat{a} ^\dagger \hat{a}$, so finally get
\begin{equation}
\hat{H} = \hbar \omega \left(\hat{N} + \frac{1}{2} \right) .
\end{equation}
Now go aside a bit and consider creation and annihilation operators. By definition,
\begin{equation}
\hat{a}^\dagger \left| n \right\rangle = \sqrt{n + 1} \left| n + 1 \right\rangle ,
\end{equation}
\begin{equation}
\hat{a} \left| n \right\rangle = \sqrt{n} \left| n - 1 \right\rangle ,
\end{equation}
where $\left| n \right\rangle$ is the eigenstate of creation and annihilation operators, as well as of the Hamiltonian (due to the fact that they commute - homework to prove).
Now
\begin{equation}
\hat{a}^\dagger\hat{a} \left| n \right\rangle = \hat{a}^\dagger \sqrt{n} \left| n - 1 \right\rangle = \sqrt{n} \sqrt{n} \left| n \right\rangle = n \left| n \right\rangle ,
\end{equation}
so conclude that the eigenvalue of a number operator, $\hat{N}$, is just $n$, so if we now apply Hamiltonian in the Schroedinger equation, get
\begin{equation}
\hat{H} \psi = E \psi ,
\end{equation}
\begin{equation}
E_n = \hbar \omega \left(n + \frac{1}{2} \right) ,
\end{equation}
which is exactly the result you were looking for.
Answer 2:
First of all you should remember that the general aim of solving an eigenvalue problem is to find a set of eigenvectors, but not a single eigenvector. In your case, equation should be modified to
\begin{equation}
\frac{d^2 \psi_n}{dx^2} + \left[(2n + 1) - \varepsilon^2\right]\psi_n = 0 ,
\end{equation}
where $\psi_n$ are eigenvectors (eigenfunctions) that correspond to eigenvalues $E_n$. Try to think a little bit and explain physical meaning of having many energy eigenvalues in quantum mechanics.
Now return to the general theory of eigenvalue equations. Although I have never met the equation you wrote, I cannot find any place it can be wrong apart from the one just pointed out. Though, I don't see how far can you go from it.
Answer 3:
Hermite polynomials are usually beyond standard quantum mechanics courses. If you know Legendre, Chebyshev and/or other polynomials, you may guess that Hermite polynomials are derived as solution to some differential equation, and this does not contradict to the definition of $\psi$.
As I've already mentioned, Hermite polynomials are usually beyond standard quantum mechanics courses. Usually you are not supposed to derive them at this level. However, if you are still interested, you may want to consult with google or ask another question here.
Hope your questions have now been answered in full. However, should you need any further comment - you are welcome.
The key concept to look for is displaced number states. These are, quite simply, the number states $|n⟩$, moved by the displacement operator
$$D(\alpha)=\exp\left[\alpha a^\dagger-\alpha^*a\right]$$
to some point $\alpha=x+ip$ on the complex phase space. The ground state of a harmonic oscillator which has been displaced to a real $\alpha=x$ is, as you note, the coherent state
$$|\alpha⟩=D(\alpha)|0⟩.$$
The rest of the eigenstates are, similarly,
$$|\alpha,n⟩=D(\alpha)|n⟩.$$
The simplest way to prove this is to use the commutation/displacement relations between the displacement operator and the ladder operators,
$$aD(\alpha)=D(\alpha)(a+\alpha)\quad\text{and}\quad a^\dagger D(\alpha)=D(\alpha)(a^\dagger+\alpha^*).$$
You can then transform the eigenvalue equation $a^\dagger a|n⟩=n|n⟩$ into
$$
\begin{align}
(a^\dagger-\alpha^*)(a-\alpha)|\alpha,n⟩
&=\left(D(\alpha)a^\dagger D(\alpha)^\dagger\right)
\left(D(\alpha)aD(\alpha)^\dagger\right) D(\alpha)|n⟩
\\ &= D(\alpha)a^\dagger a|n⟩=nD(\alpha)|n⟩
\\&=n~|\alpha,n⟩,
\end{align}
$$
or alternatively
$$
\left[a^\dagger a-(\alpha a^\dagger +\alpha^*a)\right]~|\alpha,n⟩
=(n-|\alpha|^*)~|\alpha, n⟩.
$$
Now, as all states in Hilbert space, the displaced number states can be written in the basis of number states, as $|\alpha,n⟩=\sum_{m=0}^\infty c_m(\alpha)|m⟩$. The coefficients in this expansion are simply the matrix elements of the displacement operator in the number basis:
$$⟨m|\alpha,n⟩=⟨m|D(\alpha)|n⟩.$$
This can be calculated in a couple of ways, which either look messy or seem pulled out of thin air, but what matters is ultimately the result. They come out as
$$
⟨m|D(\alpha)|n⟩=\sqrt{\frac{n!}{m!}}\alpha^{m-n}e^{-\tfrac12|\alpha|^2}L_n^{(m-n)}(|\alpha|^2)\quad\text{when }m\geq n,
\tag 1
$$
where $L_n^{(m-n)}$ is a Laguerre polynomial.
The standard reference in the literature for this matrix element is Appendix B of
Ordered Expansions in Boson Amplitude Operators. K. E. Cahill and R. J. Glauber. Phys. Rev. 177 no. 5, 1857-1881 (1969).
It's worth noting that most of the literature simply quotes the result $(1)$ and attributes it to Cahill and Glauber, but the most papers cite both the paper above and a second back-to-back paper, also by Cahill and Glauber, which does not contain the result. So be careful and show that you've read what you cite!
My bachelor's thesis,
E. Pisanty. Estados coherentes generalizados y estructura analítica del operador de aniquilación [pdf] (Generalized coherent states and the analytical structure of the annihilation operator). Tesis de licenciatura, Universidad Nacional Autónoma de México, 2011.
was on this topic and it contains an alternative calculation, though it's all in Spanish (un?)fortunately.
Best Answer
It is almost correct, but you're missing a factor of $\hbar\omega$. Note that in your equation $\partial_t b_H = \frac{\mathrm i}{\hbar} [H,b]_H$, the units do not match (since $H$ and $b$ are dimensionless). Correct is $$ \partial_t b_H = \frac{\mathrm i}{\hbar} [h,b]_H = \mathrm i\omega\, [H,b]_H . $$
Otherwise, just one hint for making your life easier: $$ [H,b]_H = [H_H,b_H] = [H, b_H] . $$