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 online. In the Appendix, "A Numerical computation of the bounce solution", it outlines an algorithm involving a shooting method to solve a boundary value problem.
The value of the constants in the equation are as follows:
Mp = 2.435*10^18;
? = 1.4*10^-5;
? = 6.3*10^-8;
? = -0.013;
The appendix outlines the algorithm, where Eq.(2.23) is referring to the differential equation (and boundary conditions) stated above.
This is my attempt thus far:
zero = 10^-10;
inf = 10^9;
ref = inf - 10^3;
\[Delta] = 0.01;
DU = Re[D[
1/4 y[x]^4*(\[Gamma] + \[Alpha]*(Log[y[x]])^2 + \[Beta]*(Log[
y[x]])^4), y[x]]];
eqn = y''[x] + 3 y'[x]/x == DU;
B = 2.3;
n = 0;
k = 0;
\[Epsilon] = 10^-10;
\[CapitalDelta] = 1;
\[CapitalLambda] = 1;
Bvalues = {};
absdiffvalues = {};
diffvalues = {};
\[Delta]values = {};
eqn0 = B + B^3/8 (? + (?/2)*Log[B] + ?*(Log[B])^2 + ?*(Log[B])^3 + ?*(Log[B])) x^2;
While[\[CapitalLambda] > \[Epsilon] && n <= 500,
eqn0 = B +
B^3/8*(\[Gamma] + (\[Alpha]/2)*
Log[B] + \[Alpha]*(Log[B])^2 + \[Beta]*(Log[
B])^3 + \[Beta]*(Log[B])^4)*x^2;
initialcon1 = eqn0 /. x -> zero;
initialcon2 = D[eqn0, x] /. x -> zero;
sol = First[
NDSolve[{y''[x] + 3 y'[x]/x == DU, y[zero] == initialcon1,
y'[zero] == initialcon2}, y[x], {x, zero, inf},
Method -> {"ExplicitRungeKutta"}]];
x2sol = x^2*y[x] /. sol;
\[CapitalDelta] = (x2sol /. x -> inf) - (x2sol /. x -> ref);
\[CapitalLambda] = Abs[\[CapitalDelta]];
AppendTo[diffvalues, \[CapitalDelta]];
AppendTo[absdiffvalues, \[CapitalLambda]];
If[n >= 3,
B0 = Part[Bvalues, n];
B1 = Part[Bvalues, n - 1];
B2 = Part[Bvalues, n - 2];
If[! (B0 > B1 > B2 || B2 > B1 > B0),
k++;
\[Delta] = \[Delta]/10;
];
];
If[\[CapitalDelta] > \[Epsilon],
B += \[Delta];
];
If[\[CapitalDelta] < -\[Epsilon],
B -= \[Delta];
];
AppendTo[Bvalues, B];
n++;
If[Divisible[n, 500],
Print[n];
Print[\[CapitalDelta]];
Print[B];
Print[k];
Print[\[Delta]];
Print[initialcon1];
Print[initialcon2];
];
];
Unfortunately it seems to get to a value of delta = -0.000129013 but then seems to not get any better no matter how many iterations it runs through. In the paper it seems as if it only takes around 100 iterations in this case to reach a suitable solution. Is there something I am missing, i.e. to search for appropriate initial guess, changing the numerical method?
Any help on this matter would be greatly appreciated.