# Solving a boundary value problem numerically, with high precision

In the paper ON A PAINLEVÉ-TYPE BOUNDARY-VALUE PROBLEM, the authors consider the BVP given by the ODE $$y''=y^2-x \tag{1}$$ with the boundary conditions \begin{align} y(0)&=0, \tag{2a} \\ y(x)& \sim \sqrt{x} \text{ as }x \to \infty \tag{2b}.\end{align} They do so by studying the 1-parametric family of solutions $y_\alpha(x)$ of (1) with the initial conditions \begin{align} y(0)&=0, \tag{3a} \\ y'(0)&=\alpha. \tag{3b} \end{align} It is conjectured* that the family $y_\alpha$ is classified into 5 different categories, depending on the value of $\alpha$ (please see attached plots):

1. If $\alpha > 0.924376$, $y_\alpha$ increases monotonically, and blows up in finite time (two "uppermost" red curves).
2. If $\alpha \approx 0.924376$, $y_\alpha$ increases monotonically, and approaches $\sqrt{x}$ from below as $x \to \infty$ (upper green curve).
3. If $-3.79199<\alpha<0.924376$, $y_\alpha$ oscillates about, and approaches $-\sqrt{x}$ as $x \to \infty$ (both blue curves).
4. If $\alpha \approx -3.79199$, $y_\alpha$ goes down first, attains a minimum, and then goes up and approaches $\sqrt{x}$ from below as $x \to \infty$ (lower green curve).
5. If $\alpha<-3.79199$, $y_\alpha$ goes down first, attains a minimum, and then intersecting $\sqrt{x}$ from below, resulting in a finite time blowup (two "lowermost" red curves).

I'm interested in finding those "critical slopes", which solve the BVP (1)-(2) (namely about 0.924376 and -3.79199). I don't understand how the authors obtained their decimal approximation to the upper critical slope, so I've decided to try and find them myself in the following way:

In order to find the upper slope for example, I have used Mathematica's numerical ODE solver. I noticed that using a slope of $1$ results in a blowup, and using a slope of $0$ results in oscillatory behaviour. From there I continued bisecting the interval $[0,1]$, closing in on the right slope. My problem with this method is that I'm at the mercy of the error tolerance of "NDSolve" and I'm not sure how good my approximation really is.

Here is my question: What is a numerical method that can find both solutions of the BVP (1)-(2) to high precision (that is, determines both $\alpha$s up to $10^{-16}$, or preferably even better) efficiently?

To be clear, I am not interested as much in the actual solution $y_\alpha(x)$, the value of $\alpha$ itself suffices.

Thank you!

$*$ Part of this classification problem has been proven in a later paper.

• did you try Taylor expanding y at x=0? Maybe one can find a recursion formula for the coefficients. Then you can restrict to this recursion instead of the ODE. Jan 15, 2017 at 21:13
• @user35593 It is proven that all solutions of this equation have a pole at some negative point. In fact there exists a number $L$, such that all real solutions have a pole at $[-L,0)$. This bounds the radius of convergence of the Taylor series. Jan 16, 2017 at 0:06