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

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!

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}]


• @W.Robin This isn't the first time I met a equation that using the asymptotic solution as the boundary can make the solving easier. (But I never encountered one that can't be solved successfully without the asymptotic solution, actually I've struggled for more than half an hour before I noticed this!) Jan 9 '16 at 8:03
• @ xzczd, 1. your answer employed NDSolve in which the 4 ICs have been estimated by power series expanding near singularity point in a consistent way. Thus, it essentially is a numerical power series method instead of the shooting method. Am I right? 2. How can I figure out the Method used by the NDSolve. I found that if I explicitly use Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}, there are some zigzag curves, so the default method is better. Thanks! Jan 9 '16 at 10:41
• @W.Robin 1. I think this method still belongs to shooting method, the power series is just to find a more accurate approximation of $y′(x_0)$ and $y‴(x_0)$. 2. This post may be helpful, BTW, if you simply set Method -> "ExplicitRungeKutta" the resulting plot will be as good as the default. Jan 10 '16 at 2:45
• @W.Robin Check this post. Jan 11 '16 at 5:45
• @W.Robin 1. FindRoots2D is essentially a function that uses points on the plot of equation as starting points to find the accurate roots of equation with FindRoot so it should be as effective as "FindRoot with a bunch of good starting points", if you still feel worried, just plug the roots back to the equation to check if they're correct . 2. Check this post. BTW, you may also give a try to other answers under the FindAllCrossings2D question. Jan 11 '16 at 15:17