RKNT86[f_List, vars_List, initialvalues_List, dinitialvalues_List,
        finalx_, errorTolerance_] :=
Block[{x = SetAccuracy[First[initialvalues], 33],
   y = SetAccuracy[Rest[initialvalues], 33], dy = SetAccuracy[dinitialvalues, 33],
   xend = SetAccuracy[finalx, 33], h, err, hnew, y7, y8, dy7, dy8, solution, hmax,
   ireject, k = Table[i + j, {i, 9}, {j, Length[vars] - 1}],
   b=SetAccuracy[{46704396222138759/1124501888012545693,0,
        84069894477030747/424535379079037893,60269691739898297/328032958547368465,
            2009963068113133/27794099874007722,162341471393132/140140455957185117,
                                    6086576956589044/1882413506280312633,0,0},33],
   bb=SetAccuracy[{10769958754260247/261191895425614637,0,
         104933541030533329/527807735255158343,8187542127950603/44863180380403502,
         50493885750265423/674323734860213804,-396215365808089/252398506959352750,
                                    5468871271464350/1319483122963052413,0,0},33],
   db=SetAccuracy[{46704396222138759/1124501888012545693,0,
       90371972523959954/390135632629351589,118990880894033457/367654647557162744,
      180830119624415039/624884373647391279,16628088200566168/1643600751401035359,
    1524820183138666476/417332398303375801,-942444174868320016/265473221553563103,
                                                                           0},33],
   dbb=SetAccuracy[{10769958754260247/261191895425614637,0,
       58861559987617091/253105545276009947,142913350550568712/444546485690175277,
          8398007711885933/28026591338889651,-8440103966850896/615634893567208211,
                                           1592393294195924241/339999309740023022,
                               -6699802037196600096/1421037300124099357,3/20},33],
   c=SetAccuracy[{0,8065253268/111157879849,16130506536/111157879849,99/229,
                                             1855/2473,116/131,1129/1130,1,1},33],
   a=SetAccuracy[{{0},{502615833312847/190946037812928939},
        {1601030787675953/456179150746555700,1601030787675953/228089575373277850},
      {47478115875661981/518814108724307373,-64883723802385428/357040639400014459,
                                            25666007926449694/139746227660637731},
    {-328112826298039228/251912779790891183,969895830706346953/297412056373654755,
    -958305119264262743/492487831928632961,151603443293999467/564549369158251216},
{44079989458325648760/345626831710945999,-267609305840442666747/859338149021870938,
 130442442641184422881/655209191357439877,-7381158156698807543/475346800759815547,
                                           594932629852457670/835908452635682287},
{-10802627635977292643/544607328597417370,22047268993379696720/454307750813938153,
   -9705881798108421635/315306127829247354,1078781161885226048/413453123878982063,
         -8616008188673363/388077019471353686,365346507915481/466435620062528214},
{-13306779498890004275/660225117657805349,22208114914951831801/450387553598953907,
   -6398475501845852180/204556450443208783,1412284034546646006/533270054097053815,
           -19179472816466775/820785347597843378,14435103384615/18331075303513484,
                                                -364401779978/904202609357507829},
    {46704396222138759/1124501888012545693,0,84069894477030747/424535379079037893,
          60269691739898297/328032958547368465,2009963068113133/27794099874007722,
  162341471393132/140140455957185117,6086576956589044/1882413506280312633,0}},33],
   half = SetAccuracy[1/2, 33], two = SetAccuracy[2, 33],
   ninetenth = SetAccuracy[9/10, 33],
   ode = Function[Release[vars], Release[f]]},
   hmax = SetAccuracy[xend - x, 33]; h = SetAccuracy[errorTolerance^(1/8), 33];
   ireject = 0;
   solution = SetAccuracy[{Flatten[{x, y, dy}]}, 33];
   k[[9]] = SetAccuracy[Apply[ode, Flatten[{x, y}]], 33];
   While[x < xend, If[x + h == x, Break[]];
      If[x + h - xend > 0, h = SetAccuracy[xend - x, 33]];
      k[[1]] = k[[9]];
      k[[2]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[2]] h, y + c[[2]] h*dy
                                              + h^2*a[[2, 1]] k[[1]]}]], 33];
      k[[3]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[3]] h, y + c[[3]] h*dy
                                             + h^2*a[[3]].Take[k, 2]}]], 33];
      k[[4]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[4]] h, y + c[[4]] h*dy
                                             + h^2*a[[4]].Take[k, 3]}]], 33];
      k[[5]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[5]] h, y + c[[5]] h*dy
                                             + h^2*a[[5]].Take[k, 4]}]], 33];
      k[[6]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[6]] h, y + c[[6]] h*dy
                                             + h^2*a[[6]].Take[k, 5]}]], 33];
      k[[7]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[7]] h, y + c[[7]] h*dy
                                             + h^2*a[[7]].Take[k, 6]}]], 33];
      k[[8]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[8]] h, y + c[[8]] h*dy
                                             + h^2*a[[8]].Take[k, 7]}]], 33];
      k[[9]] = SetAccuracy[ Apply[ode, Flatten[{x + c[[9]] h, y + c[[9]] h*dy
                                             + h^2*a[[9]].Take[k, 8]}]], 33];
      y8 = SetAccuracy[y + h*dy + h^2 b.k, 33];
      y7 = SetAccuracy[y + h*dy + h^2 bb.k, 33];
      dy8 = SetAccuracy[dy + h db.k, 33];
      dy7 = SetAccuracy[dy + h dbb.k, 33];
      err = SetAccuracy[Max[{Max[Abs[y8 - y7]/10], Max[Abs[dy8 - dy7]/10]}], 33];
      hnew = SetAccuracy[ Min[SetAccuracy[hmax, 33],SetAccuracy[ h/Max[half, Min[two,
             SetAccuracy[(Rationalize[err/errorTolerance, 10^-33]^(1/7))/ninetenth,
             33]]], 33]], 33];
      If[err <= errorTolerance, x += h; y = y8; dy = dy8;
            AppendTo[solution, SetAccuracy[Flatten[{x, y8, dy8}], 33]],
                                  hnew = SetAccuracy[Min[hnew, h], 33];
            ireject = ireject + 1;
      k[[9]] = k[[1]]];h = hnew]; Return[{solution,ireject}]]