By the simple form of the equation (1) you wrote down, Allan really meant
$$ \frac{\partial}{\partial t} \Psi(x,t)|_{t=0} = \frac{E}{i\hbar}\Psi(x,0) $$
He just used the notation where $t=0$ is substituted from the beginning but he clearly did mean that $\Psi(x)$ is first considered as a general function of $t$, then differentiated, and then we substitute $t=0$.
This equation says that the time derivative of $\Psi(x,t)$ at $t=0$ is proportional to the same wave function. By itself, it does not imply that $\Psi(x,t)$ for an arbitrary later $t$ will be given by equation (2): if we only constrain the derivative at one moment $t=0$, the wave function may do whatever it wants at later (or earlier) moments $t$.
However, we may generalize (1) to any moment $t$ which is what you wrote down
$$ \frac{\partial}{\partial t} \Psi(x,t) = \frac{E}{i\hbar}\Psi(x,t) $$
and this equation does imply (2). If the $t$-derivative of $\Psi(x,t)$ is proportional to the same $\Psi(x,t)$, then $\Psi(x,t)$ and $\Psi(x,t')$ are proportional to each other for each $t,t'$. That implies that $\Psi(x,t)$ must factorize to $\Psi(x)f(t)$ and the function $t$ is the simple complex exponential to solve the equation with the right coefficient.
Because you are more or less writing the same things, I find it plausible that you are not missing anything.
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.
Best Answer
Those substitutions are quite helpful. Not because of the addition of new information (as you mentioned, we don't learn anything which we didn't know before), but because of the way they let us see what we already have. They also let you eyeball the system, as I like to call it, which essentially means that you can guess a lot more without resorting to computational tools.
Essentially, after making these substitutions,
You can easily think about what would happen if you change the proportions between the input values (this is a general advantage of nondimensionalization)
You can 'guess' possible solutions easily
Here's an explanation:
You didn't go through the substitutions, so here's a brief overview of the process of solving it. You can check each stage to make sure that you see what exactly is happening to the dimensions of the equation. Start off with the Schrödinger equation, $$\left(\frac{P^2}{2m}+\frac{1}{2}\omega^2x^2\right)\psi(x, t)=E\psi(x, t)$$ (Note that the LHS is the Hamiltonian, $\mathscr{H}=T+V=\left(\frac{P^2}{2m}+\frac{1}{2}\omega^2x^2\right)$, and this is all clearly dimensionally correct)
You're familiar with the concept of non-dimensionalizing, so I'll skim over that: take $\epsilon=\frac{E}{\hbar\omega}$, and get (after expanding your momentum operator and stuff) $$\frac{m\omega}{2\hbar}x^2\psi(x, t)-\frac{\hbar}{2m\omega}\frac{\partial^2}{\partial x^2}\psi(x, t)=\epsilon\psi(x, t)$$ Right now, there's no obvious difference between this equation and the original with $E$ in it: although it looks a bit cleaner, there's nothing trivially obvious when you look at it. So let's try another pair of substitutions, $x=\alpha u$ and $\alpha=\sqrt{\frac{\hbar}{m\omega}}$. Most books show these substitutions made one after the other, with expansions of the equation at each stage, but since the question says that they're already mathematically self-explanatory, I'll skip through to the result of applying those substitutions and simplifying stuff. You get this magically clean expression, $$u^2\psi(x, t)-\frac{\partial^2}{\partial u^2}\psi(x, t)=2\epsilon\psi(x, t)$$ Or, if I can be a bit looser with the notation, $$-\psi''+u^2\psi=2\epsilon\psi$$ and $$\psi''=(u^2-2\epsilon)\psi$$
This is obviously quite easy to solve in certain conditions. We can guess what happens if $u\rightarrow\infty$: the $\epsilon$-related terms become negligibly small, so we solve $\psi''=u^2\psi$. And that's easy enough to guess; the results are along the lines of $\psi=Au^ke^{u^2/2}$. Clearly that's not something you could achieve if you didn't have $u$ in the picture, and if you hadn't introduced $\epsilon$, it would have been sticky to figure out exactly what terms were being approximated to zero. A similar process can be followed to approximate solutions for $u\rightarrow 0$.
The next obvious advantage of removing the units is cleanness and generality of computations for different orders of magnitude. It's a common technique to assume all other units in a system to take values such that $\hbar=1$ and $c=1$. It makes a lot of calculations easier. In this case, we're setting the units of energy ($E$) as $\hbar\omega$ (just to verify, $\hbar$ has units of action, $[E\ T]$ (I'm using a nonstandard $[E]$ for energy), and $\omega$ is frequency, $[T^{-1}]$). We're also scaling the length to $\sqrt{\frac{\hbar}{m\omega}}$. Additionally, by setting these values, it's easy to make intuitive guesses about the scale of the problem (this is a more general advantage of removing dimensions, and isn't very specific to the quantum harmonic oscillator). This makes it easy to use the same equation for systems with different orders of magnitude (for example situations where the amplitude $u$ is really tiny and situations where it's huge), without loosing an understanding and 'feeling' for what's going on.