Finding solutions of the differential equation $x\frac{d}{d x}\left(x\frac{d y}{d x}\right)=4(y-x^{2}y+y^{3})$ In my research I have come across the following non-linear differential equation:
$$x\frac{d}{d x}\left(x\frac{d y}{d x}\right)=4(y-x^{2}y+y^{3})$$
I want to find the general solution of this equation for $x\geq0$, but didn't manage to do it.
I have found some facts about the solutions by investigating the asymptotics:

*

*For $x\ll1$, the solutions go like $y=Cx^2$ for some constant $C$, or diverge.

*For $x\gg1$, the solutions go like $y=x$, or go to $0$ (probably after oscillating, found by numerics), or diverge.

In particular, I am interested in a one of the solutions. By looking upon the numerics, I believe that there is only one solution that goes like $y=x$ for $x\gg1$, and it does not diverge at $x=0$. I didn't manage to prove this by now.
I am interested in finding this solution analytically, or if impossible, just find appropriate $C$ analytically.
Thanks!
 A: Mathematica's command DSolve cannot solve your ODE in closed form; so, apparently, a closed-form solution does not exist.
However, it will shown here that, in line with your conjecture, for each real $c$ your ODE has solutions $y$ such that $y(x)=(c+o(1))x^2$ as $x\to0+$.
Indeed, your ODE is
\begin{equation*}
    x^2 y''+x y'=4(y-x^2y+y^3) \tag{1}
\end{equation*}
for real $x>0$.
Assuming the initial conditions $y(0+)=y'(0+)=0$ and using the Taylor formulas $y'(x)=x\int_0^1 v(sx)\,ds$  and $y(x)=x^2\int_0^1 (1-s)v(sx)\,ds$ with $v:=y''$, we rewrite (1) as
\begin{equation*}
    v(x)=\int_0^1 K(x,s)v(sx)\,ds+4x^4\Big(\int_0^1 (1-s)v(sx)\,ds\Big)^3, \tag{2}
\end{equation*}
where
$$K(x,s):=[4(1-x^2)(1-s)-1]s.$$
For any real $c$, setting now
\begin{equation*}
    v(x)=c+x w(x),
\end{equation*}
we rewrite (2) as
\begin{equation*}
    w(x)=F(w)(x):=-2cx+\int_0^1 K(x,s)w(sx)\,ds+4x^3\Big(\int_0^1 (1-s)(c+xw(sx))\,ds\Big)^3. \tag{3}
\end{equation*}
Take now any real $m>0$ and $h>0$, and let
\begin{equation*}
    W_{m,h}:=\{w\in C[0,h]\colon\|w\|\le m\},
\end{equation*}
where $\|w\|:=\max_{x\in[0,h]}|w(x)|$.
Since $|K(x,s)|$ is convex in $x^2$, for $x\in[0,1]$ we have
\begin{equation*}
    \int_0^1 |K(x,s)|\,ds\le\max\Big(\int_0^1 |K(0,s)|\,ds,\int_0^1 |K(1,s)|\,ds\Big)
    =\max(19/48,1/2)=1/2. 
\end{equation*}
So, for any $w\in W_{m,h}$
\begin{equation*}
\|F(w)\|\le2|c|h+m/2+h^3(|c|+hm)^3/2\le m \tag{4}
\end{equation*}
if $h>0$ is small enough -- which will be assumed in what follows. So, $F$ maps $W_{m,h}$ into $W_{m,h}$.
Moreover, if for some $w$ and $u$ in $W_{m,h}$ we have $\|w-u\|\le t$ for some real $t\ge0$, then similarly to (4) we get
\begin{equation*}
    \|F(w)-F(u)\|\le t/2+3h^3(|c|+hm)^2ht/2\le\tfrac23\,t
\end{equation*}
if $h>0$ is small enough. So, $F$ is a contractive map of $W_{m,h}$ into $W_{m,h}$, and hence $F$ has a fixed point. That is, equation (3) has a solution $w\in W_{m,h}$.
Thus, ODE (1) has a solution $y$ with $y''(x)=c+x w(x)$ for $x\in[0,h]$ and $w\in W_{m,h}$, and with $y(0+)=y'(0+)=0$. So, $y(x)=cx^2/2+O(x^3)$ as $x\to0+$.
$\Box$
A: Is this too naive? To obtain the large-$x$ behavior, I take a series expansion of $y(x)$ in powers of $1/x$,
$$y(x)=c x + a_0 + a_1/x + a_2/x^2+ a_3/x^3 +\cdots,$$
substitute into
$$R(x)=x\frac{d}{d x}\left(x\frac{d y}{d x}\right)-4(y-x^{2}y+y^{3}),$$
expand $R(x)$ in powers of $1/x$ and demand that the leading order terms vanish. It follows that $c=1$, $a_{n}=0$ for any even $n$, while for odd $n$ I find
$$a_1=-\frac{3}{8},\;\;a_3=-\frac{9}{128},a_5=-\frac{99}{1024},\;\;a_7=-\frac{11637}{32768},\;\;a_9=-\frac{627669}{262144},$$
so coefficients of the form $a_{2n+1}=-b_n/2^{2n+1}$ with integer $b_n$.
