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):

- If $\alpha > 0.924376$, $y_\alpha$ increases monotonically, and blows up in finite time (two "uppermost" red curves).
- If $\alpha \approx 0.924376$, $y_\alpha$ increases monotonically, and approaches $\sqrt{x}$ from below as $x \to \infty$ (upper green curve).
- If $-3.79199<\alpha<0.924376$, $y_\alpha$ oscillates about, and approaches $-\sqrt{x}$ as $x \to \infty$ (both blue curves).
- 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).
- 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.