Storing results of a calculation as a list

Clear[Pressure, EnergyDensity, fromPressureToEnergy, 
  fromEnergyToPressure, Ae, Be, Ap, Bp, Er, Mr, Pr, r, h, g1, f1, k1, 
  m1, f2, g2, k2, m2, f3, g3, k3, m3, f4, m4, g4, k4, M, R]; 

Ae = Import["D:Energy density.xlsx", {"Data", 1, All, 1}];
Be = Ae*(10^-6/1.1154907);
EnergyDensity = Be/(0.08969);
Ap = Import["D:Pressure.xlsx", {"Data", 1, All, 1}];
Bp = Ap*(10^-6/1.1154907);
Pressure = Bp/(0.08969);
fromPressureToEnergy = 
  Interpolation[Most@Transpose[{Pressure, EnergyDensity}], 
   InterpolationOrder -> 1];
fromEnergyToPressure = 
  Interpolation[Most@Transpose[{EnergyDensity, Pressure}], 
   InterpolationOrder -> 1];
Er = (1637.99731*(10^-6/1.1154907)/0.08969); (*Initial value*)
h = 0.0001;
Mr = 10^-33;
Pr = fromEnergyToPressure[Er];
r = h;

While[Pr > 0, 

 g1 = 4*[Pi]*r^2*Er*.08969;
 f1 = -(1.47*Mr*
      Er*(1 + Pr/Er)*(1 + 4*[Pi]*r^3*Pr*.08969/Mr))/(r^2*(1 - 
       2*Mr*1.47/r));
 k1 = h*f1;
 m1 = h*g1;
 
 f2 = -(1.47*Mr*
      Er*(1 + (Pr + k1/2)/Er)*(1 + 
        4*[Pi]*(r + h/2)^3*(Pr + k1/2)*.08969/Mr))/((r + h/2)^2*(1 - 
       2*Mr*1.47/(r + h/2)));
 g2 = 4*[Pi]*(r + h/2)^2*(Er + m1/2)*.08969;
 k2 = h*f2;
 m2 = h*g2;
 
 f3 = -(1.47*Mr*
      Er*(1 + (Pr + k2/2)/Er)*(1 + 
        4*[Pi]*(r + h/2)^3*(Pr + k2/2)*.08969/Mr))/((r + h/2)^2*(1 - 
       2*Mr*1.47/(r + h/2)));
 g3 = 4*[Pi]*(r + h/2)^2*(Er + m2/2)*.08969;
 k3 = h*f3;
 m3 = h*g3;
 
 f4 = -(1.47*Mr*
      Er*(1 + (Pr + k3)/Er)*(1 + 
        4*[Pi]*(r + h)^3*(Pr + k3)*.08969/Mr))/((r + h)^2*(1 - 
       2*Mr*1.47/(r + h)));
 g4 = 4*[Pi]*(r + h)^2*(Er + m3)*.08969;
 k4 = h*f4;
 m4 = h*g4;

 Pr = Pr + 1/6*(k1 + 2*k2 + 2*k3 + k4);
 Mr = Mr + 1/6*(m1 + 2*m2 + 2*m3 + m4);

 Er = fromPressureToEnergy[Pr];
 r = r + h;
 ]

M = Mr; (*result*)
R = r;(*result*)

Dear Friends in the code above I have an initial value for (Er) which is (1637.99731*(10^-6/1.1154907)/0.08969) and the final results are (M) and (R).

My question is that how I can solve this code for various (Er) which varies from the initial value above to (135*(10^-6/1.1154907)/0.08969) and for each of the (Er) I want to save M and R. Therefore, finally I could have something like a list that consists of different (M) and (R) that related to various (Er) which are considered as the initial value.

In fact, my code only works for one (Er) and as a result one (M) and (R) for the chosen (Er). What I need to do is solving my code for various (Er) regarding the mentioned interval.