ClearAll[DirectSmoothCubicFit];

DirectSmoothCubicFit[y_, m_] :=
 Module[{n, F1, F2, F, a, b, k, , exp1, exp2, exp3, exp4, exp5, exp6,
   eqns, c, coffarray, coff, finalvals = {}},
  n = Length[y];
  F1 = Sum[(y[[j]] - a[0] - a[1]*(j - m) -
        Sum[a[k]*(j - m)^k, {k, 2, 3}])^2, {j, 1, m}]/2;      (*
  Definition of the Error Function *)
  F2 = Sum[(y[[j]] - a[0] - a[1]*(j - m) -
        Sum[b[k]*(j - m)^k, {k, 2, 3}])^2, {j, m, n}]/2;
  F = F1 + F2;

  exp1 = D[F, a[0]];   (* [PartialD]F/[PartialD]a[0] *)
  exp2 = D[F, a[1]];   (* [PartialD]F/[PartialD]a[1] *)
  exp3 = D[F, a[2]];   (* [PartialD]F/[PartialD]a[2] *)
  exp4 = D[F, a[3]];   (* [PartialD]F/[PartialD]a[3] *)
  exp5 = D[F, b[2]];   (* [PartialD]F/[PartialD]b[2] *)
  exp6 = D[F, b[3]];   (* [PartialD]F/[PartialD]b[3] *)

  eqns = {exp1 == 0, exp2 == 0, exp3 == 0, exp4 == 0, exp5 == 0,
    exp6 == 0};
  {c, coffarray} =
   Normal[
    CoefficientArrays[
     eqns, {a[0], a[1], a[2], a[3], b[2], b[3]}]];    (*
  Calculation of the Coefficient Matrix for the six equations *)

  coff = LinearSolve[coffarray, -c];
  finalvals =
   Join[finalvals,
    Table[coff[[1]] + coff[[2]]*(j - m) + coff[[3]]*(j - m)^2 +
      coff[[4]]*(j - m)^3, {j, 1, m}]];
  finalvals =
   Join[finalvals,
    Table[coff[[1]] + coff[[2]]*(j - m) + coff[[5]]*(j - m)^2 +
      coff[[6]]*(j - m)^3, {j, m + 1, n}]];

  finalvals
  ]

then the loop where the problem is shown

y =  RandomVariate[NormalDistribution[0, 1]]*
    Table[i^3, {i, 1, 20}] + RandomVariate[NormalDistribution[0, 1]];
fact = 1./0.000001^2;
For[o = 5, o <= 16, o++, fit1 = DirectSmoothCubicFit[y, o];
 residuals1 = Abs[(fit1)] // N ;
 res1 = Total [residuals1^2]*fact // N;
 r1 = Join[r1, {res1}]];
r1