numerical integration - Solving a second-order nonlinear differential equation

I am trying to solve a particular Cauchy problem given by

I found from a particular paper that the solutions looks like

For the auxiliary conditions

For only specific values of $a_{i}$

I found from another paper that the problem can be solved with the Mathematica code

Off[General::stop]

createLargestRegion[h_, d_, mult_] :=


 Block[{ bound = Min[mult*Abs[d]/h, 19], n = 1, result, previous }, 
  While[Not[
     Check[result[r_] := 
       Evaluate[
        First[\[Theta][r] /.

          NDSolve[{\[Theta]''[r] + (1/r ) \[Theta]'[r] -
              (0.5/(r^2)) Sin[
                2 \[Theta][r]] == ((h /2.0) Sin[\[Theta][r]] + 
               (d /r) Sin[\[Theta][r]]^2 ), \[Theta][0.001] == 

             Pi, \[Theta][bound] == 0.1}, \[Theta][r], r, 
           "Method" -> {"Shooting", 
             "StartingInitialConditions" -> { \[Theta][0.001] == 
                Pi, \[Theta]'[0.001] == (-Pi/10)}}]]], False]] || 
    FindMaximum[{result[t], t > 0.1 && t <= bound}, t] > Pi || 
    result [0.5] > result [0.001], bound = bound - 0.1;
   n++;];
  result[r]]

hdInterpFuncs = 

 Table[Table[{i*0.1, j*-0.1, 
    Evaluate[createLargestRegion[i*0.1, j*-0.1, 19.0]]}, {j, 3, 18, 
    3}], {i, 3, 18, 3}]

I completeley understand NDSolve and Shooting method, but i am unable to understand the rest of the code. Also when i run the code i am getting

NDSolve::berr: The scaled boundary value residual error of 6.892078251140068`*^7 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

Any help in helping me understand the code and or figuring out why i am getting the error will be thankful. Also if there are other ways to solve the problem with better code will be thankful too.

(English is not my first language forgive grammar errors.)

Related Paper here