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.