Trouble with shooting method for a 4th-order differential equation

Question

I'm trying to solve the following forth-order ODE with the shooting method:

$$\frac{1}{5}(y-2xy^\prime)=\frac{1}{x}\left\{\frac{xy^\prime}{y}+xy^3 \left[\frac{(xy^\prime)^\prime}{x} \right]^\prime \right\}^\prime$$ where $\prime$ denotes differentiation. It is clear that this ODE need 4 boundary conditions (BCs): two of which are given $y^\prime(0)=y^{\prime\prime\prime}(0)=0$, however, the other two values of $y(0)$ and $y^{\prime\prime}(0)$ are determined by shooting for BCs at infinity ($x_\text{Max}$).

It has been found that $y=Ax^{1/2}$ is a far-field asymptotic solution to the ODE by neglecting the nonlinear terms. Now, let me refer to the solutions with $y\sim x^{1/2}$ asymptotic behavior as the $x^{1/2}$ solutions. Here, I want to search just for this $x^{1/2}$ solutions by shooting method in which the ODE is integrated with a 4th-order R-K scheme from the origin to a certain end-value of $x_{max}$. In practice, I truncate the domain to $x_0<x<x_\text{Max}$ to avoid the singularity at $x=0$, and I would like to work with $x_0=10^{-4}$ and $x_\text{Max}=10$.

The shooting parameters $y(0)$ and $y^{\prime\prime}(0)$ will be adjusted so that $y-2xy^\prime=0$ or $y+4x^2y^{\prime\prime}=0$ at $x=x_\text{Max}$. The two BCs were chosen to be independent to involve low-order derivatives and to require $y\propto x^{1/2}$ at the end of the integration interval.

My questions are:

(1) How can I use a Taylor series (say, five-terms) to start the integration at $x_0=10^{-4}$. I am thinking the series should be located in the position of initial conditions, but I can't figure it out;

(2) How do I implement this searching method: firstly, find the solutions that satisfying the 1st BCs $y-2xy^\prime=0$ at $x_\text{Max}$, next, find the solutions that satisfying the 2nd BC $y+4x^2y^{\prime\prime}=0$ at $x_\text{Max}$. Then ParametricPlot the shooting parameters $(y(0), y^{\prime\prime}(0))$ that corresponding to solutions satisfying the above-mentioned BCs respectively, thus the intersections of the two sets of curves correspond to the desired $x^{1/2}$ solutions.

(3) Is it possible to use Do loop and FindRoot to find the parameters, as shown in this answer?

Thanks!

Answer 1

Interesting equation. It seems to be necessary to use the asymptotic solution as the boundary at $x=x_0$ if you want to solve the equation correctly.

Thanks to this answer, we can easily get the series solution at $x=0$ with seriesDSolve:

eq1 = 1/5 (y[x] - 2 x y'[x]) == D[(x y'[x])/y[x] + x y[x]^3 D[D[x y'[x], x]/x, x], x]/x; 
bc1 = y[x] - 2 x y'[x] == 0; 
bc2 = y[x] + 4 x^2 y''[x] == 0; 

solSeries = seriesDSolve[eq1, y, {x, 0, 5}, 
                         {y[0] -> a, y'[0] -> 0, y''[0] -> b, y'''[0] -> 0}]    

With the asymptotic solution, we obtain a list of more accurate boundaries at $x=x_0$:

newbclist = Thread[(Derivative[#][y][x0] == 
        (D[Normal@solSeries, {x, #}] /. x -> x0) & ) /@ Range[0, 3]]

Finally, solve the equation and plug it into the right boundary and plot:

x0 = 1/10^4; xMax = 10;
sol = ParametricNDSolveValue[{eq1, newbclist}, y, {x, x0, xMax}, {a, b}]

ContourPlot[{bc1, bc2} /. x -> xMax /. y -> sol[a, b] // Evaluate, 
            {a, 0.4, 0.8}, {b, -0.05, 0.25}]

By contrast, if one simply approximate the left boundary as $y'(x_0)=y'''(x_0)=0$, the resulting plot will be undesired:

bclist = Thread[(Derivative[#][y][x0] & /@ Range[0, 3]) == {a, 0, b, 0}]

solfake = ParametricNDSolveValue[{eq1, bclist}, y, {x, x0, xMax}, {a, b}]

ContourPlot[{bc1, bc2} /. x -> xMax /. y -> solfake[a, b] // Evaluate, 
            {a, 0.4, 0.8}, {b, -0.05, 0.25}]