arf2[rpk_, jpka_, n_Integer] := Module[{rk1, rk2, s}, 
     {rk1, rk2} = rpk; Return[arf1b[f, {rk1, 1, rk2}, jpka, n]]; ]
 
arf1b[f_, regp_, jpk_, m_Integer] := Module[{r, s, syn, aa, aacc, jw, kw, 
      cht}, {jw, kw} = jpk; cht = {zz21, jw, kw}; s = 0; Label[o1]; 
      {r, cht} = rconti[f, regp, cht]; s = s + 1; If[s < m, Goto[o1]]; 
      {syn, aa, aacc} = r; Return[syn]; ]
  
rconti[f_, regpk_, chpk_] := Module[{sy1a, a, acc, jv, kv, r, ch1, ch2, ch3, 
      s}, {sy1a, a, acc} = regpk; {ch1, jv, kv} = chpk; 
      ch2 = rf1c[f, sy1a, a, acc]; ch3 = aver[ch1, ch2, jv, kv]; 
      r = {sy1a + ch3, a, acc}; s = {ch3, jv, kv}; Return[{r, s}]; ]
 
rf1c[f_, sy_, a_, n_Integer] := (AccuracyGoal -> n; PrecisionGoal -> n; 
     WorkingPrecision -> n + 7; Module[{u, zy, du, b, fe, fe2, st21, st22}, 
      b = a; st22 = N[sby[{{D[f, y1], D[f, y2], D[f, y3], D[f, y4], D[f, y5], 
            D[f, y6], D[f, y7], D[f, y8], D[f, y9], D[f, y10], D[f, y11], 
            D[f, y12], D[f, y13], D[f, y14], D[f, y15], D[f, y16], D[f, y17], 
            D[f, y18], D[f, y19], D[f, y20], D[f, y21]}, f}, sy], n]; 
       st22 = Rationalize[st22, 1/10^(n + 2)]; {st21, fe} = st22; 
       du = sqsum[st21]; u = -(fe/du); Goto[o2]; Label[o1]; b = (2*b)/3; 
       Label[o2]; zy = N[sy + b*u*st21, n]; fe2 = N[sby[f, zy], n]; 
       If[fe2 < fe, Goto[o3]]; If[b < 9^(-n), Return[gonetoosmall]]; 
       Goto[o1]; Label[o3]; Return[zy - sy]; ])
  
aver[chy1_, chy2_, j_Integer, k_] := Module[{v}, v = j*chy1 + chy2; 
      v = v/(1 + j); v = k*v; Return[v]; ]