Optimizing a slow For loop

I know that this is a "beginner question", and, in fact, I am. I want to improve my code, as it takes really too much time to run.
I have already read some other discussions, like here,
but I am struggling in translating the simplifications with Table or Do into my example.

The For cycle I want to improve is the following:

zh = 0.4;
list = {};

For[i = 1, i < 10000, i++, 
    With[{zg = -0.6, dh = i*10^-2}, 
         nsol = Block[{eps = $MachineEpsilon}, 
         NDSolve[{phi''[x] + 2*phi'[x]/x + (2*(zg + zh*Exp[-zh*x/dh])/x + 1)*phi[x] == 0, phi[eps] == 1, phi'[eps] == -(zg + zh)}, phi, {x, eps, 20000}, 
                 WorkingPrecision->MachinePrecision, AccuracyGoal->15, PrecisionGoal->8, MaxSteps->Infinity]]];  
    AppendTo[list, 1/Evaluate[(15000*phi[15000])^2 + ((15000-Pi/2)*phi[15000-Pi/2])^2 /. nsol[[1]]]];]

Clearly, this code, written in this way, is highly inefficient. Also, I need to do more of these, with different values for zg inside With, and make some plots out of the lists.

Anyone that can help me with this noob question?
Thanks a lot!

phifunc = ParametricNDSolveValue[
  {[Phi]''[x] + 2 [Phi]'[x]/x + (2 (zg + zh Exp[-zh x/dh])/x + 1) [Phi][x] == 0
   , [Phi][$MachineEpsilon] == 1
   , [Phi]'[$MachineEpsilon] == -(zg + zh)
   }
  , 1/((15000 [Phi][15000])^2 + ((15000 - [Pi]/2) [Phi][
         15000 - [Pi]/2])^2)
  , {x, $MachineEpsilon, 20000}
  , {dh [Element] Reals, zg [Element] Reals, zh [Element] Reals}
]
list = ParallelTable[phifunc[i 10^-2, -0.6, 0.4], {i, 1, 10000}]