I am looking to replicate an algorithm explained in the paper "Impact of new physics on the EW vacuum stability in a curved spacetime background" by E. Bentivegna, V. Branchina, F. Continoa and D. Zappal` to numerically solve the boundary-value problem encountered when trying to find the bounce solution in false vacuum decay. It can be found here: https://arxiv.org/abs/1708.01138
In the Appendix, "A Numerical computation of the bounce solution", it outlines an algorithm involving a shooting method to solve the following boundary problem:
$y^{\prime\prime}(x) + \frac{3}{x}y^{\prime}(x) = \frac{dU}{dy};\\ y(\infty) = 0,\ y^{\prime}(0) = 0, \\U(y) = \frac{1}{4}y^{4}(\gamma +\alpha\ln^{2}y + \beta\ln^{4}y),$
where ${}^{\prime}$ indicates differentiation with respect to $x$.
The value of the constants in the equation are as follows:
Mp = 2.435*10^18;
\[Alpha] = 1.4*10^-5;
\[Beta] = 6.3*10^-8;
\[Gamma] = -0.013;
\[Lambda]6 = 0;
The appendix outlines the algorithm as follows, where Eq.(2.23) is referring to the differential equation (and boundary conditions) stated above. 
I am attempting to recreate this algorithm in Mathematica but my attempts thus far of simply trying numerically integrate the differential equation using a shooting method with boundary conditions specified in the Appendix has just yielded errors owing to finding complex infinities.
Any help with this matter would be greatly appreciated.
Attachments:
