Although this is a question pertaining to Physics, since this is related to the mathematical treatment of a differential equation, I believe it is well suited for this community.
While deriving the wave function for harmonic oscillator potential using Schrodinger's equation, we obtain the following equation through rearrangement of constants and nondimensionalization of the variables. $$\frac{d^2\psi}{du^2}+(\epsilon-u^2)\psi=0 \tag1$$ And then we use the technique of asymptotic analysis.
This is achieved by checking the behaviour of $\psi$ at large $u$ and guessing the form of the solution as $$\psi \approx \exp(-u^2) g(u)$$
And then we obtain Hermite's differential equation for $g(u)$ which can solved by power series solution.
My Question: Why can't we avoid Asymptotic Analysis and directly go for a Series Solution? Why can't we just directly take $$\psi(u)=\sum_\limits{n=0}^{\infty} a_n u^n \tag2$$
I have checked everywhere in the internet and also in all standard books on quantum mechanics. What I have observed is that they directly go for the asymptotic analysis without stating any reason. They simply say that the asymptotic analysis will help in simplifying the calculations. However they do not mention anything about a direct solution by power series method. They do not make any comments on the possible of a direct series solution; neither why we may be able to go for such solutions nor why we can't go for such a method and have to adopt something called asymptotic analysis.
I tried solving Schrödinger's Equation using such a power series as in $(2)$. What I got was: $$2a_2+\epsilon a_0+(6a_3+\epsilon a_1)u+\sum_\limits{n=0}^{\infty} \left[(n+4)(n+3)a_{n+4}+\epsilon a_{n+2}-a_n\right]u^{n+2}=0$$
These gives 2 constants and $1$ recursion. $$(n+4)(n+3)a_{n+4}+\epsilon a_{n+2}-a_n=0$$
I know that it is difficult to obtain any nice desired result from this recursive relation.
But is this correct? Is this process feasible here? I checked for singularities and found none.
In case this is correct, is it so that both asymptotic analysis and my procedure are allowed but the asymptotic analysis method is most favoured since in that case we get closed form results which can be used to derive other useful results?
Or is this series solution not feasible due to some reason more general?
Your solution $\psi$ as wave function should be bounded, even go to zero at infinity. It should at least be in $L^2(\Bbb R)$, for the equation to make sense you will need a Sobolov space $H^1$ or better. That is not given for polynomial solutions, and also most power series solutions will diverge at infinity, trigonometric and similar functions are the exceptions to that rule.
Thus it makes sense to get the solution as the product of a factor that captures this asymptotic behavior and a slower growing factor that can be computed via power series, or in this case can be chosen as polynomial so that the convergence considerations occur in the linear combinations of the basis solutions.
You will find that in setting $\psi(u)=\exp(s(u))g(u)$ you get \begin{align} \psi''&=(e^s)''g+2(e^s)'g'+e^sg''\\ &=e^s\Bigl[(s'^2+s'')g+2s'g'+g\Bigr]\\ &=e^s(u^2-E)g \end{align} Assuming $s'$ is growing to infinity, that is, $s$ is superlinear, the terms with $s'^2$ and $u^2$ are the largest on both sides. Setting them equal (this is a design decision, there is one functional relation with $s$ and $g$ free) to cancel each other in the equation results in $s'=\pm u$, $s=\pm u^2/2$ where the plus sign gets excluded because of unbounded growth. $\psi(u)=\exp(-u^2/2)g(u)$ has the remaining equation $$ 0=(E-1)g-2ug'(u)+g''(u) $$ or in power series coefficients $$ 0=(E-1)a_n-2na_n+(n+1)(n+1)a_{n+1}\implies a_{n+2}=\frac{2n+1-E}{(n+2)(n+1)} $$ To prevent a solution growing like $e^{u^2/2}$ you want either a polynomial solution or that $\frac{a_{n+2}}{a_n}=o(\frac1n)$. The latter is impossible. For a polynomial solution you need $E=2N-1$ to get $a_n=0$ for $n>N$.