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.
Indeed, as suggested by phase-space quantization, most of these equations are reducible to generalized Laguerre's, the cousins of Hermite. As universally customary, I absorb $\hbar$, M and ω into r,E. Note your E is twice the energy.
Since $r\geq 0$ you don't lose negative values, and you may may redefine $r^2\equiv x$, so that
$$
r\partial_r = 2x \partial_x \qquad \Longrightarrow r\partial_r (r\partial_r)= r^2\partial_r^2+ r\partial_r=4(x^2\partial_x^2+x\partial_x),
$$
hence your radial equation reduces to
$$
\left ( \partial_x^2+ \frac{1}{x}\partial_x +\frac{E-x}{4x} -\frac{m^2}{4x^2} \right ) R(m,E)=0 ~.
$$
Now, further define
$$
R(m,E)\equiv x^{|m|/2} e^{-x/2} ~ \rho(m,E),
$$
to get
$$
\partial_x R(m,E)= x^{|m|/2} e^{-x/2} \left (-1/2 +\frac{|m|}{2x} + \partial_x \right )~ \rho(m,E) \\
\partial_x^2 R(m,E)= x^{|m|/2} e^{-x/2} \left (-1/2 +\frac{|m|}{2x} + \partial_x \right )^2~ \rho(m,E),
$$
whence the generalized Laguerre equation for non-negative m=|m|,
$$
x \partial_x^2\rho(m,E) +\left({m+1} -x\right )\partial_x \rho(m,E)+\frac{1}{2}(E/2-m-1) \rho(m,E)=0~.
$$
This equation has well-behaved solutions for non-negative integer $$k=(E/2-m-1)/2\geq 0 ~,$$ to wit, generalized Laguerre (Sonine) polynomials $L^{(m)}_k (x)=x^{-m}(\partial_x -1)^k x^{k+m}/k!$.
Plugging into the factorized solution and the above substitutions nets your eigen-wavefunctions. The ground state is $k=0=m$, ($E=2$ in your conventions), so a radially symmetric Gaussian, $e^{-r^2/2}$.
Again, in your idiosyncratic convention, the degeneracy is E/2.
So, degeneracy 2 for $E=4$ : $m=1$, $k=0$; you may check this is just $r e^{-r^2/2 +i\phi} $. You may choose the $\cos \phi$ and $\sin \phi$ solutions, if you wish, constituting a doublet of the underlying degeneracy group SU(2).
Best Answer
I assume you have no qualms with the "large $\xi$" approximation - it's fairly obvious that $\xi^2-k^2\approx \xi^2$ for large enough $\xi$. After that you're left with the differential equation $$\frac{d^2\psi}{d\xi^2}\approx-\xi^2\psi.\tag1$$ One way to solve this equation is by the method of divine inspiration: you somehow come up with two linearly independent functions you can write down which solve the equation, after which you know what the general solution is.
However, the two functions that Griffiths poses are not exact solutions of that equation, as you would know if you had done your proper diligences. Indeed, $$ \frac{d^2}{d\xi^2}\left[e^{\pm \xi^2/2}\right] =\frac{d}{d\xi}\left[\pm\xi e^{\pm \xi^2/2}\right] =(\xi^2\pm1)e^{\pm\xi^2/2}. $$ This means that we're looking for a solution that's approximately valid for the (already approximate) equation (1). As for how one might get such solutions, there's obviously a myriad different possible paths.
One really nice way of deriving the solutions is to factor the offending second-order differential operator into two different first-order operators: $$ \left(\frac{d}{d\xi}-\xi\right)\left(\frac{d}{d\xi}+\xi\right) \approx \left(\frac{d}{d\xi}+\xi\right)\left(\frac{d}{d\xi}-\xi\right) \approx \frac{d^2}{d\xi^2}-\xi^2. $$ These are of course approximate inequalities, and of course you must work these out to see what terms got dropped and why.
After that, you can simply work out solutions to the two equations $\left(\frac{d}{d\xi}-\xi\right)\psi=0$ and $\left(\frac{d}{d\xi}+\xi\right)\psi=0$, which are in fact the solutions given by Griffiths. These are first-order equations and therefore simply solvable by integrating. The (approximate) equalities above guarantee that a function in the kernel of either factor will be (approximately) in the kernel of the operator you do care about.
Of course, these factors are far from trivial inventions, and they are at the heart of the operator approach to solving the simple harmonic oscillator.