(*RKN 8(6) with 9 stages,no-FSAL*)
RKNT86q9[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[{3191538187421696/76607108605432915,0,13815874303602012/69579866183121917,14604812893174087/79378705834398872,12061218770183621/166622303733231213,15609617015400/233291059437933767,371765604219257/111475530824146994,0,0},33],
bb=SetAccuracy[{4544292102832777/109056534231464193,0,4682651711005479/23585400043481548,46722285954615265/253893219962912894,4751354290135738/65721585748949841,20872833551830/134159415686285343,275420922524446/83046920983443867,0,0},33],
db=SetAccuracy[{3191538187421696/76607108605432915,0,10308242332317290/44357423208271919,7107618457535881/21873268413857328,22056521909108756/75044404292647497,15596425292979/34434009875005756,325257858967320448/9895379989758637,-(264730262449877449/7963593493382224),17208373/35885750},33],
dbb=SetAccuracy[{4544292102832777/109056534231464193,0,18333976229602070/78901367072948263,146694624662575579/451359699798674378,40221502534828457/137021353651599420,91894267481143/87253900673082639,776789986225611057/23764274461164518,-(1116801360586595899/33934531992244452),23651021/71771500},33],
c=SetAccuracy[{0,2595146787461113/35654960162808999,23785164771277655/163393282122478121,14427641/33259908,26914142/35708683,15577224/18277247,38090011/38093876,1,1},33],
a=SetAccuracy[{{0},{295132092736843/111419829353054663},{378512699615967/107173587955359337,802015671331405/113542950051902326},{9945580188014483/107861941766479192,-(21127832523454066/115356389813386625),18088716445271473/97760613913942175},
               {-(184569114806220359/112841400437628580),595308873796066195/146500969503370446,-(95938071830688501/39190010187048758),24740235889975229/81328315902644410},
               {828692824853675681/365166841077510,-(8922830626242929564/1616145596072733),5309688443105545745/1512691814917754,-(1024584250889564737/3835297140363491),243334944688840544/26685249097802661},
               {-(198499310481410068/14988189920044743),988020934248343439/30631779844455146,-(372584950612767755/18396862167620476),93706617067436735/54962117018052057,2652169291282213/72706638769934851,3326767107636/45583415053986647}, 
               {-(172476446800076249/77764528330584470),1052320941122775251/32896321613528843,-(287682559714467205/6581569888910478),336649615658501777/14242861273902858,-(94884627/9749078),-(177655963/35046632),112476592/20068355},
               {1589642054066860483/2111418052567415,206513499/21728459,-(5009179395035143313/3047562608623994),2192653675860564860/1440780190602451,-(1099957025566422337/1624301323788501),-(3640940497065881569/10360892974776789),1917284830561677115/4934686172719308,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];
While[x<xend,If[x+h==x,Break[]];
If[x+h-xend>0,h=SetAccuracy[xend-x,33]];
k[[1]]=SetAccuracy[Apply[ode,Flatten[{x,y}]],33];
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]/100000],Max[Abs[dy8-dy7]/100000]}],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];
h=hnew];Return[{solution,ireject}]]




(*--------- example ------------*)
(* Solve Inho Linear system from
   from: Mathematics/MDPI, 2023, 11, 891. 
 https://www.mdpi.com/2227-7390/11/4/891
*)
{solution, ireject} = 
  RKNT86q9[{1/100*z1 - 1/10*z2, -1/10*z1 + 1/100*z2 + Sin[x]}, {x, z1,
     z2}, {0, 1, 1}, {-1000/10101, -(10100/10101)}, 10*Pi, 10^-22];