How to get solutions in elementary functions to the following non-linear ODE $y\frac{\text{d}^2y}{\text{d}x^2}-2\left(\frac{\text{d}y}{\text{d}x}\right)^2+xy^3\frac{\text{d}y}{\text{d}x}=0$. I want to get solution in elementary functions in order to do asymptotic analysis. I tried the problem with Mathematica. It did not give any solutions.
 A: HINT : 
You have to do asymptotic analysis without solving the ODE. 
If you had the analytic solution of the ODE and then study the asymptotic behavior you distort the problem.
I suppose that the ODE was especially chosen so that one cannot express the exact solution on a closed form.
If the initial conditions are given (for example at $x=0\quad\to\quad y(0)=y_0$  and $y'(0)=y'_0$ ) the function $y(x)$ is determined. Since $y(x)=C$=constant is a solution $y(x\to\infty)=C$. The function tends asymptotically to a constant (not always the same, depending on the initial conditions).
This is confirmed by numerical calculus, for examples :
 
$$y\frac{\text{d}^2y}{\text{d}x^2}-2\left(\frac{\text{d}y}{\text{d}x}\right)^2+xy^3\frac{\text{d}y}{\text{d}x}=0$$
As $x\to\infty \qquad y\to C.\quad$ So, the ODE is then approximately :
$$C\frac{\text{d}^2y}{\text{d}x^2}-2\left(\frac{\text{d}y}{\text{d}x}\right)^2+C^3 x\frac{\text{d}y}{\text{d}x}\simeq 0$$
$$C\frac{\text{d}^2y}{\text{d}x^2}
+\left(-2\frac{\text{d}y}{\text{d}x}+C^3 x \right)\frac{\text{d}y}{\text{d}x}
\simeq 0$$
Moreover, $\frac{\text{d}y}{\text{d}x}$ is small while $x$ is large. Thus :
$$\frac{\text{d}^2y}{\text{d}x^2}+C^2 x\frac{\text{d}y}{\text{d}x}
\simeq 0$$
$$y(x)\simeq C-c_1\left(1-\text{erf}\left(\frac{C}{\sqrt{2}}x\right) \right)\qquad x \text{ large}$$
A: It might be also useful to transform this equation into a system of first-order ODEs. That might be helpful for some numerical approximations and analysis of the dynamics. 
Let $y' = \frac{dy}{dx}$. Then equation is 
$$y\, y'' - 2 \, (y')^2 + xy^3 \, y' = 0$$ Multiply both sides of the equation by $\frac{1}{y^3}$ and get
$$\frac{y''}{y^2} - 2 \,\frac{(y')^2}{y^3 } + x \, y' = 0$$ Observe that 
$$\frac{d}{dx}\left(\frac{y'}{y^2}\right) = \frac{d}{dx}\left(y^{-2}\, {y'}\right) = {y^{-2}}{y''} + (-2) \, y^{-3} \,( y')^2 = \frac{y''}{y^2} - 2 \,\frac{(y')^2}{y^3 }$$ leading to the equation
$$\frac{d}{dx}\left(\frac{y'}{y^2}\right) + x\, y' = 0$$ Now, set 
$$z = \frac{y'}{y^2}  \,\,\,\,\ y' = y^2 \, z$$ This gives you the first-order system
\begin{align}
 \frac{dy}{dx} &= y^2  z\\
\frac{dz}{dx} &= - x \, \frac{dy}{dx} = -x \, y^2 z
\end{align}
Finally it is a non-autonomous first order system of ODEs 
\begin{align}
 \frac{dy}{dx} &= y^2  z\\
\frac{dz}{dx} &= -x \, y^2 z
\end{align}
If you write it as an autonomous system in the extended phase space, by introducing $x(t) = x_0 + t\, ,$ then
\begin{align}
 \frac{dx}{dt} &= 1\\
 \frac{dy}{dt} &= y^2  z\\
\frac{dz}{dt} &= -x \, y^2 z\\
\end{align}
which is equivalent to the vector field 
$$X = \frac{\partial}{\partial x} +  y^2  z \, \frac{\partial}{\partial y}   -x \, y^2 z \, \frac{\partial}{\partial z} $$
Observe that $X$ is tangent to the field of two dimensional planes in three space defined by the contact one form
$$\omega = x \, dy + dz$$ Also notice that the coordinate planes $y=0$ and $z=0$ are invariant planes and the solutions on these are constants, i.e. $x(t)= x_0 + t, \, y(t) = y_0, \, z(t) = 0$ and  $x(t)= x_0 + t, \, y(t) = 0, \, z(t) = z_0$.
Maybe this is helpful to run some analysis to figure out the behavior of the solutions.
Edit. Here is how the equation is derived as well as some related equations. We are looking for a fucntion $\theta(s,t)$ satisfying the diffusion equation $$\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial s}\left(D(\theta) \, \frac{\partial \theta}{\partial s}\right)$$ for some function $D(\theta)$ and for some boundary value data. Then one looks for solution of the form $$\theta(s, t) = y\big(s \, t^{-1/2}\big)$$ for some function $y=y(x)$. Then, having in mind that $x = s \, t^{-1/2}$, 
$$\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial t} \,  y\big(s \, t^{-1/2}\big) = -\,\frac{1}{2}\,\frac{dy}{dx}\,s\, t^{-3/2} =  -\,\frac{1}{2}\,s\, t^{-3/2} \,\frac{dy}{dx}$$
$$\frac{\partial \theta}{\partial s} = \frac{\partial}{\partial s} \,  y\big(s \, t^{-1/2}\big) = \frac{dy}{dx}\, t^{-1/2}= t^{-1/2} \, \frac{dy}{dx}$$ so the equation becomes
$$\frac{\partial \theta}{\partial t} = -\,\frac{1}{2}\,s\, t^{-3/2} \,\frac{dy}{dx} = \frac{\partial}{\partial s} \left(D(y) \, t^{-1/2} \, \frac{dy}{dx}\right)$$
$$-\,\frac{1}{2}\,s\, t^{-3/2} \,\frac{dy}{dx} = \frac{\partial}{\partial s} \left(D(y) \, t^{-1/2} \, \frac{dy}{dx}\right) = D'(y)\, t^{-1} \left(\frac{dy}{dx}\right)^2 + D(y)\, t^{-1} \, \frac{d^2y}{dx^2}$$ so dividing by $t^{-1}$ both sides leads to 
$$-\,\frac{1}{2}\,s\, t^{-1/2} \,\frac{dy}{dx} = D'(y)\, \left(\frac{dy}{dx}\right)^2 + D(y) \, \frac{d^2y}{dx^2}$$
Recall that $x = s \, t^{-1/2}$ so
$$-\,\frac{1}{2} \, x \,\frac{dy}{dx} = D'(y)\, \left(\frac{dy}{dx}\right)^2 + D(y) \, \frac{d^2y}{dx^2}$$ so the equation becomes
$$D'(y)\, \left(\frac{dy}{dx}\right)^2 + \, D(y) \, \frac{d^2y}{dx^2} \, + \,\frac{1}{2} \, x \,\frac{dy}{dx} = 0$$ and can be rewritten as
$$\frac{d}{dx}\left(  \, D(y) \, \frac{dy}{dx} \, \right) + \,\frac{1}{2} \, x \,\frac{dy}{dx} = 0$$ or if you prefer 
$$\frac{d}{dx}\left(  \, 2 \, D(y) \, \frac{dy}{dx} \, \right) + \, x \,\frac{dy}{dx} = 0$$ 
Your equation arises when $$D(y) = \frac{1}{2 \, y^2}$$
Then one can set $z = 2\, D(y) \frac{dy}{dx} $, so $ \frac{dy}{dx} = \frac{z}{2\,D(y)} $ leading to the system of ODEs
\begin{align}
\frac{dy}{dx} &= \frac{z}{2\, D(y)} \\
\frac{dz}{dx} &=  -\, x \,\frac{dy}{dx} = -\,  \frac{x \, z}{2 \,D(y)}
\end{align} or in autonomous form
\begin{align}
\frac{dx}{d\tau} &= 1\\ 
\frac{dy}{d\tau} &= \frac{z}{2 \, D(y)} \\
\frac{dz}{d\tau} &= -\, \frac{x \, z}{2\,D(y)}
\end{align} 
There some other manipulations one can carry out like
$\frac{d}{dx}(xy) = x \frac{dy}{dx} + y$ so the equation is also 
$$\frac{d}{dx}\left( \, 2 \, D(y) \, \frac{dy}{dx}  + x\,y \, \right) \, - \, y = 0$$ Furthermore, one can also interpret $x$ as a function of $y$, that is $x = x(y)$ so $$\frac{d}{dx}\left(  \, 2 \, D(y) \, \frac{dy}{dx} \, \right) +  x \,\frac{dy}{dx} = 0$$ can be written as
$$\frac{dy}{dx}\frac{d}{dy}\left( \, 2 \, D(y) \, \frac{dy}{dx} \, \right) \, + \, x \,\frac{dy}{dx} = 0$$ so
$$\frac{d}{dy}\left( \, 2 \, D(y) \, \frac{dy}{dx} \, \right) \, +  \, x  = 0$$
$$\frac{d}{dy}\left(  \, \frac{ 2\, D(y)}{\frac{dx}{dy}} \, \right) \, + \, x  = 0.$$ Now, if you set $\varphi(y)$ such that $\varphi'(y) = x(y)$ then
$$\frac{d}{dy}\left( \, 2 \, \frac{ D(y)}{\varphi''} \, \right) \, +  \, \varphi'  = 0$$ and if you integrate once, up to constant
$$\frac{2\, D(y)}{\varphi''} \, +  \, \varphi  = 0$$ equivalently
$$D(y) = \,\frac{1}{2} \, \varphi \, \frac{d^2\varphi}{dy^2}$$
The solvability mentioned in the article you referenced, comes from simply taking appropriate function $\varphi(y)$ and plugging them in the latter equation to obtain expressions for $D(y)$. That's how the author has created his table there. Studying the symmetries of any of the mentioned avatars of the same equation is a much harder task. Probably more interesting. I unfortunately do not see an obvious symmetry in the case $D(y)= \frac{1}{2y^2}$, which is the case of the equation is your posting. I am not sure about the case $D(y) = \sqrt{y}$ either.   
Edit. Reduction and idea for explicit solution. I found a way to solve the equation. I guessed a one parameter Lie symmetry whose vector field had first integral and I used that first integral to simplify the coordinates, so that the equation becomes solvable. 
Consider the one parameter (parameter $\sigma \in \mathbb{R}$) Lie group acting on the space $(x,y)$
\begin{align}
x &\mapsto e^{-\sigma} x\\ 
y &\mapsto e^{\sigma} y  
\end{align} 
 This group leaves the equation 
\begin{align}
\frac{d}{dx}\left(\frac{1}{y^2} \, \frac{dy}{dx} + x y \right) - y = 0
\end{align} 
invariant, so it is a symmetry, mapping a solution of the equation to a solution of the eqution. The vector field that generates this Lie gourp symmetry is simply
$$Z_L = -\,x\,\frac{\partial}{\partial x} + y\,\frac{\partial}{\partial y}$$
which has a first integral, i.e. conserved quantity, invariant quantity, etc...
$$J(x,y) = xy$$
This tells us that one can use $J$ in order to find an explicit change of variables that will transform the vector field $Z_L$ into the straight, constant vector field $\frac{\partial}{\partial \tau}$. This means that in the new coordinates, one of the variables in our second order differential equation will disappear (it becomes a cyclic coordinate, i.e. explicitly non-existent). Here is what the coordinates should be
\begin{align}
\tau &= - \log(x)\\
u &= J(x,y) = xy\\
\end{align}
Then the inverse change of variables is
\begin{align}
x &= e^{-\tau}\\
y &= e^{\tau} u
\end{align} 
The differentials of this change of coordinates are
\begin{align}
dx &= - \, e^{-\tau} \, d\tau\\
dy &= e^{\tau} u \, d\tau + e^{\tau} \, du
\end{align} 
We plug in the equation in its form
\begin{align}
\frac{d}{dx}\left(\frac{1}{y^2} \, \frac{dy}{dx} + x y \right) - y = 0
\end{align} 
and we obtain
\begin{align}
- \,\frac{d}{e^{-\tau}d\tau}\left(\frac{1}{e^{2\tau}u^2} \, \frac{e^{\tau} u \, d\tau + e^{\tau} \, du}{e^{-\tau}d\tau}  + u\right) - e^{\tau}u = 0
\end{align} 
All the exponents $e^{\tau}$ miraculously cancel out
\begin{align}
\frac{d}{d\tau}\left(\frac{1}{u^2} \, \frac{u \, d\tau +  du}{d\tau}  + u\right) + u = 0
\end{align} 
leaving us with an autonomous equation!
\begin{align}
\frac{d}{d\tau}\left(\frac{1}{u}  + \frac{1}{u^2} \,  \frac{du}{d\tau}  + u\right) + u = 0
\end{align} Moreover, this equation can be simplified further. Make another change of variables
$$w = \frac{1}{u}$$ Then 
$$\frac{dw}{d\tau} = - \frac{1}{u^2} \, \frac{du}{d\tau}$$
so the latter differential equation transforms into
\begin{align}
\frac{d}{d\tau}\left(w  - \frac{dw}{d\tau}  + \frac{1}{w}\right) + \frac{1}{w} = 0
\end{align}
Introduce the extra variable
$$v = w  - \frac{dw}{d\tau}  + \frac{1}{w}$$ which leads us
the equivalent system of first order ODEs
\begin{align}
\frac{dw}{d\tau} &= w   + \frac{1}{w} - v\\
\frac{dv}{d\tau} &= - \frac{1}{w}
\end{align}
After eliminating $d\tau$ this becomes
$$\frac{dw}{w +\frac{1}{w} - v } = w \, dv$$ or equivalently
$$\big(w^2 +1 - v \, w\big)  \, dv = dw$$
Observe that we have arrived at a first order ordinary differential equation. If we solve it, we will have a formula that links $v$ to $w$ and thus a first order autonomous ODE for $w$ only, which is solvable. Indeed, the latter ODE is 
$$\frac{dw}{dv} = w^2 +1 - v \, w$$ 
$$\frac{dw}{dv}  - 1 = w^2 - v \, w$$
$$\frac{d}{dv} \big(w-v\big)  = w\,(w -v)$$ If we set $W = w-v$ then $w = W + v$ and thus
$$\frac{dW}{dv} = v\,W + W^2$$ which is a standard solvable equation. I think it is called Bernoulli type equation and is linearizible.  
Indeed your original equation is solvable. However, I am going to stop here and continue later.  
This is a very nice illustration of how Lie group theory in the form of symmetries lead to the reduction of the problem to a simpler one. 
