Here is my updated code to see if this helps people. It still doesn't yield the required results but the algorithm follows the same procedure as in the paper.
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 = 0.070;
n = 0;
k = 1;
\[Epsilon] = 10^-10;
\[CapitalDelta] = 1;
\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\) = 1;
Bvalues = {};
absdiffvalues = {};
While[
\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\) > \[Epsilon],
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 = NDSolve[{y''[x] + 3 y'[x]/x == DU, y[inf] == 0,
y'[zero] == 0}, y[x], {x, zero, inf},
Method -> {"BoundaryValues" -> {"Shooting",
"StartingInitialConditions" -> {y[zero] == initialcon1,
y'[zero] == initialcon2}}}];
x2sol = x^2*y[x] /. (First@sol);
\[CapitalDelta] = (x2sol /. x -> inf) - (x2sol /. x -> ref);
\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\) =
Abs[\[CapitalDelta]];
AppendTo[absdiffvalues,
\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\)];
If[n >= 3,
B0 = Part[Bvalues, n];
B1 = Part[Bvalues, n - 1];
B2 = Part[Bvalues, n - 2]
];
If[\[CapitalDelta] > \[Epsilon],
If[B1 <= B0 <= B2 ,
k++;
B += \[Delta]/(k*10);
];,
B += \[Delta]/(k*10);
];
If[\[CapitalDelta] < -\[Epsilon],
If[B2 <= B0 <= B1 ,
k++;
B -= \[Delta]/(k*10);
];,
B -= \[Delta]/(k*10);
];
AppendTo[Bvalues, B];
n++;
If[Divisible[n, 500],
Print[n];
Print[\[CapitalDelta]];
Print[B];
Print[k];
Print[initialcon1];
Print[initialcon2];
];
];
Any help with this problem would still be greatly appreciated.
