Hello guys, I've been working on ParametricNDSolveValue with 4th order Runge-Kutta method and sometimes I received this error

ParametricNDSolveValue::mxst: Maximum number of 10000 steps reached at the point x$214882$215402 == 0.9256478123350554`.

InterpolatingFunction::dmval: "Input value {1.} lies outside the range of data in the interpolating function. Extrapolation will be used."

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances

For the RK4 method my StepSize->1/1000, I've tried changing the PrecisionGoal and AccuracyGoal, only the point of x changes but also with the same error. This error does not occur all the time, if the I change the value of the inputs, sometimes I don't receive this error

Below is the code from other threads I had post and received help from the community.

prog := Module[{Twater = 300., Tnew = 0, Cwater, Dwater, Ewater, 
    ClassicalRungeKuttaCoefficients, f, x, y, point, ans, tol = 0.001,
     iters = 0, iterLim = 100, prevTnew = 10^5, 
    plotList = {}}, $MinPrecison = 50;
   Cwater = 
    QuantityMagnitude[
     ThermodynamicData["Water", 
      "ThermalConductivity", {"Temperature" -> 
        Quantity[Twater, "Kelvins"]}]];
   ans = "Twinkle, twinkle, little star,
             How I wonder what you are.
             Up above the world so high,
             Like a diamond in the sky.
             Twinkle, twinkle, little star,
             How I wonder what you are!";
   While[(Abs[prevTnew - Tnew] > tol) && (iters++ < iterLim), 
    Dwater = (5 Cwater)/Pi;
    Ewater = Dwater^3 + 2*Cwater;
    ClassicalRungeKuttaCoefficients[4, prec_] := 
     With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
       bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
      N[{amat, bvec, cvec}, prec]];
    f = ParametricNDSolveValue[{SetPrecision[
        Derivative[1][y][x] == 
         Piecewise[{{(y[x] + x^3 + 3 z - 120*Ewater), 
            0 <= x <= 1}, {(y[x] + x^2 + 2 z), 
            1 <= x <= 2}, {(y[x] + x + z), 2 <= x <= 3}}], 40], 
       y[0] == 0}, y[3.], {x, 0., 3.}, z, 
      Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
        "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
      StartingStepSize -> 1/1000];
    point = {z /. 
       FindRoot[f[z] == 100., {z, 0.1}, Evaluated -> False], 100.};
    AppendTo[plotList, point];
    Print[point];
    FindRoot[f[z] == 100., {z, 0.1}, Evaluated -> False];
    prevTnew = Tnew;
    Tnew = 
     170 + 1.89*z /. 
      FindRoot[f[z] == 100., {z, 1}, Evaluated -> False];
    Cwater = 
     QuantityMagnitude[
      ThermodynamicData["Water", 
       "ThermalConductivity", {"Temperature" -> 
         Quantity[Tnew, "Kelvins"]}]];(*Note Tnew here!*)
    Print["Twater: ", Twater, " Tnew: ", Tnew];];(*End while*)
   ans = {"Tnew: ", Tnew, "Cwater: ", Cwater, "Dwater: ", Dwater, 
     "Ewater: ", Ewater};
   Print["The answer is:\n", ans];];

The question is, I have find the StepSize to 1/1000, so wouldn't ParametricNDSolveValue run the command with step size of 0.001? Why would it stop before that? NOTE: this is the example of the model, I cannot post my real model.