(* Mathematica package that generates a 9-stages, FSAL, 8th order RKKN method *)
RKN8[g4_, g5_, g6_, g7_, a82_, db9_] :=
 Module[{e, b, db, db1, db3, db4, db5, db6, db7, db8, g, g2, g3, a21, a31, a32,
          a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72,
          a73, a74, a75, a76, a81, a83, a84, a85, a86, a87, a, mond, simp1, gg,
                                               ii, equs, so, ae, ag, ag2, bag},
e = {1, 1, 1, 1, 1, 1, 1, 1, 1};
db = {db1, 0, db3, db4, db5, db6, db7, db8, db9};
g = {0, g2, g3, g4, g5, g6, g7, 1, 1};
gg = DiagonalMatrix[g]; ii = IdentityMatrix[9]; b = db.(ii - gg);
a = {{0, 0, 0, 0, 0, 0, 0, 0, 0},
     {a21, 0, 0, 0, 0, 0, 0, 0, 0},
     {a31, a32, 0, 0, 0, 0, 0, 0, 0},
     {a41, a42, a43, 0, 0, 0, 0, 0, 0},
     {a51, a52, a53, a54, 0, 0, 0, 0, 0},
     {a61, a62, a63, a64, a65, 0, 0, 0, 0},
     {a71, a72, a73, a74, a75, a76, 0, 0, 0},
     {a81, a82, a83, a84, a85, a86, a87, 0, 0}, b};
mond = {db.e == 1, db.g == 1/2, db.g^2 == 1/3, db.g^3 == 1/4, db.g^4 == 1/5,
                                                     db.g^5 == 1/6, db.g^6 == 1/7};
simp1 = {(b.a)[[2]] == 0, (db.a)[[2]] == 0, (db.(g^2*a))[[2]] == 0};
equs = {db.(gg - ii).(gg - g7*ii).a.(gg - g3*ii).(gg - g4*ii).g ==
                Integrate[(x - 1)*(x - g7)*Integrate[Integrate[(x - g3)*(x - g4)*x,
                                               {x, 0, x}], {x, 0, x}], {x, 0,  1}],
        db.(gg - ii).a.(gg - g3*ii).(gg - g4*ii).(gg - g5*ii).g ==
                Integrate[(x - 1)*Integrate[Integrate[(x - g3)*(x - g4)*(x - g5)*x,
                                                {x, 0, x}], {x, 0, x}], {x, 0, 1}],
        db.(gg - ii).a.(gg - g3*ii).(gg - g4*ii).(gg - g6*ii).g ==
                Integrate[(x - 1)*Integrate[Integrate[(x - g3)*(x - g4)*(x - g6)*x,
                                               {x, 0, x}], {x, 0, x}], {x, 0, 1}]};
g3 = (15 - 20*g4 - 20*g5 + 28*g4*g5 - 20*g6 + 28*g4*g6 + 28*g5*g6 - 42*g4*g5*g6
          - 20*g7 + 28*g4*g7 + 28*g5*g7 -      42*g4*g5*g7 + 28*g6*g7 - 42*g4*g6*g7
                  - 42*g5*g6*g7 + 70*g4*g5*g6*g7)/(2*(10 - 14*g4 - 14*g5 + 21*g4*g5
          - 14*g6 + 21*g4*g6 + 21*g5*g6 - 35*g4*g5*g6 - 14*g7 + 21*g4*g7 + 21*g5*g7
            -35*g4*g5*g7 + 21*g6*g7 - 35*g4*g6*g7 - 35*g5*g6*g7 + 70*g4*g5*g6*g7));
g2 = g3/2;
so = Solve[mond, {db1, db3, db4, db5, db6, db7, db8}];
{db1, db3, db4, db5, db6, db7, db8} = so[[1, All, 2]];
ae = a.e - g^2/2; ag = a.g - g^3/6; ag2 = a.g^2 - g^4/12; a32 = g3^3/6/g2;
so = Solve[{ag[[4]] == 0, ag2[[4]] == 0}, {a42, a43}]; {a42, a43} = so[[1, All, 2]];
so = Solve[simp1, {a72, a62, a52}]; {a72, a62, a52} = so[[1, All, 2]];
so = Solve[equs, {a65, a76, a75}];  {a65, a76, a75} = so[[1, All, 2]];
so = Solve[Join[ag[[5 ;; 7]], ag2[[5 ;; 7]]] == Array[0 &, 6]
                                                  , {a53, a54, a63, a64, a73, a74}];
{a53, a54, a63, a64, a73, a74} = so[[1, All, 2]];
bag = db.(a - gg^2/2 + gg - ii/2);
so = Solve[bag[[3 ;; 7]] == Array[0 &, 5], {a83, a84, a85, a86, a87}];
{a83, a84, a85, a86, a87} = so[[1, All, 2]];
so = Solve[ae[[2 ;; 8]] == Array[0 &, 7], {a21, a31, a41, a51, a61, a71, a81}];
{a21, a31, a41, a51, a61, a71, a81} = so[[1, All, 2]];
Return[{a, g, b, db}]]


(*-------------------------------------------------------------------------------------------*)
(*----------------------- coefficients of RKN8(5)3 triplet ----------------------------------*)
a={{0,0,0,0,0,0,0,0,0},
   {4552000/1779408839,0,0,0,0,0,0,0,0},
   {4786441/1403290107,5099029/747467372,0,0,0,0,0,0,0},
   {154355229/1747859738,-211394569/1207216589,1079924893/6088903920,0,0,0,0,0,0},
   {-1909252503/1589450314,1717805592/570806317,-314365424/175661333,232924059/932301775,0,0,0,0,0},
   {21946941630/1188347359,-51402534543/1142361320,34478280149/1193453181,-2004368142/947175379,110757191/866705166,0,0,0,0},
   {-23246929975/1521737838,49742027226/1331069443,-10244269133/433259647,6333125575/3106146311,-305982/499127563,7096255/1058685153,0,0,0},
   {-6251522306/342681487,47239/1059,-58765771569/2079513332,1839671917/766441435,-6869719/563046564,4308559/595743943,-81644/562506817,0,0},
   {26770331/653556116,0,481915997/2463777720,83699121/458874562,83897329/1176822726,4694793/705442414,2924104/945955131,0,0}};
g={0,63781483/891695059,95796413/669639401,509/1196,669/913,1019/1180,793/798,1,1};
b={26770331/653556116, 0, 481915997/2463777720, 83699121/458874562, 83897329/1176822726, 4694793/705442414, 2924104/945955131, 0, 0};
bb={6092579/188948739, 0, 34025767/175712602, 205729963/782269670, -65978983/254934332, 47615623/122054613, 674071/220266358, -99/800, 1/2000};
bbb={33595973/765550025,0,192976793/996552243,98217979/543912306,28567777/404766130,8455799/1283407508,2867149/936899033,0,1/600};
db={26770331/653556116,0,358154341/1569106712,372703631/1173713952,223504194/837852631,45155994/925772753,238514129/483458217,-299632595/786663961,-19/1288};
dbb={13011125/280406193,0,64889674/303007291,66407981/195830378,142069414/620381939,21546492/286681337,50245305/82044331,-51725033/101893770,-2855444648/331872649939};
dbbb={39739233/916775494,0,97497477/431459905,224243619/713318590,133441962/492846883,5047045/97769056,482080661/997391952,-198326731/525952178,-7393685/599938819};

(* bb stands for 5th order, bbb stands for 3rd order, db stands for b', dbb for \hat{b}', etc *)
(*--------------------------------------------------------------------------------------------*)