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.
Best Answer
Answer to part 1: This is a common method of solving differential equations, employed by physicists to quickly extract a solution without having to make slow progress using more rigorous methods. The first step is to look for an asymptotic solution (i.e. the limit of large $\xi$ in this case), where the equation is easily solvable. Now, we know that this will not give us an exact solution, so we can't just proceed to claim $$\psi (\xi)=Ae^{-\xi^2/2}$$ However, since this solution must be exact in the limit, we can try to look for a solution of the form $$\psi(\xi)=\phi(\xi)e^{-\xi^2/2}$$ Which, we hope, may allow us to recast the original equation in a simpler or well-known form. Really, we've done nothing but use a change in variables, motivated by the asymptotic solution. Ideally, one would like to rigorously justify things like this (why exactly is it reasonable to believe that this makes things easier?) but, being a good physicist, we just push on and hope for the best ;-)
Answer to part 2: This is quite simple. To see that this is true, just try to take the limit $j\to \infty$. Since 1, 2 and $K$ are all constants, these eventually become negligible, so that we may write $$a_{j+2}=\frac{2j}{j^2}a_j$$ Which, after canceling a $j$, gives your formula. To solve this for $a_j$, follow ACuriousMind's advice (watch out: You cannot apply the recursion formula all the way to $a_0$, because $2/0$ is undefined) and realize that \begin{align*}a_3&=2a_1,\ a_4=a_2\\ a_5&=\frac{2}{3}a_3=\frac{2}{3}*2a_1,\hspace{1cm}a_6=\frac{1}{2}a_2\\ a_7&=\frac{2}{5}*\frac{2}{3}*\frac{2}{1}a_1\hspace{1cm}a_8=\frac{2}{6}*\frac{2}{4}*\frac{2}{2}a_2\end{align*}
$$ a_{2j+1}=\frac{2^j}{(2j+1)(2j-1)\dots1}a_1\hspace{1cm}a_{2j}=\frac{2^{j-1}}{2j(2j-2)\dots 2}a_2 $$ Now, for even terms $a_{2j}$, this is quite easily solved:$$a_{2j}=\frac{2^{j-1}}{2^jj!}a_2=\frac{a_2}{2}\frac{1}{j!}\implies a_k=\frac{C}{(k/2)!},\hspace{1cm}k=2j$$ Here, I used the fact that, for even $k=2n$, the double factorial $k!!$ can be easily expressed in terms of $n!$: $$k!!=2n(2n-2)*\dots *2=2(n)*2(n-1)*\dots*2(1)=2^n n!$$
For the odd terms, everything gets a bit messier. We first use the odd-numbered ($k=2n+1$) expression for $k!!$ in terms of $n!$: $$k!!=(2n+1)(2n-1)*\dots * 1=\frac{(2n+1)(2n)(2n-1)*\dots*1}{(2n)(2n-2)*\dots*2}=\frac{(2n+1)!}{2^n n!} $$ Plugging in this result, we have: $$a_{2j+1}=\frac{2^{2j}j!}{(2j+1)(2j)!}a_1 $$ Regrettably, it seems that this expression is not easy to evaluate (in terms of elementary functions). You'll need to use the gamma function instead, and most intuition is lost... Regardless, I hope that my treatment was sufficient to give you some intuition.