%                   coefsRKN.m
% Coefficients of selected RKN methods
%  S_h = (a_{s+1} b_s a_s ... b_1 a_1)
% Methods with number of stages s={1 4,6,6,7,11,14}
% There are non-symmetric compositions as well as symmetric ABA and BAB 
%  compositions: sa=1 calls for a symmetric ABA scheme and sa=2 for the BAB
% S. Blanes, F. Casas and A. Murua, Splitting methods for differential equations,
%   Acta Numerica (2024)     (Section 8, Table 8.3)

function [a b] = coefsRKN(s,sa)
% s = Number of stages 
if (s==1)  % 2nd-order: Strang
     a=[1/2 1/2];
     b=[1 0];
elseif (s==4)  % 4-stage non-symmetric 4th-order: McLachla-Atela (MA92)

    a1 = 0.5153528374311229364;
    a2 = -0.085782019412973646 ;
    a3 = 0.4415830236164665242 ;
    a4 = 1-(a1+a2+a3); %=0.1288461583653841854 ;

    b1 = 0.1344961992774310892 ;
    b2 =-0.2248198030794208058 ;
    b3 = 0.7563200005156682911 ;
    b4 = 1-(b1+b2+b3); %= 0.3340036032863214255;
    a=[0  a1 a2 a3 a4];
    b=[b1 b2 b3 b4  0];

elseif (s==6) 
    if sa==1  % 6-stage symmetric 4th-order BAB: Blanes&Moan (BM02)
	b1= 8.298440641740515e-2;
	b2= 0.396309801498368;
	b3= -3.905630492234859e-2 ;
	b4=1-2*(b1+b2+b3);
	a1= 0.245298957184271;
	a2=0.604872665711080;
	a3=1/2-(a1+a2);
    a=[0  a1 a2 a3 a3 a2 a1];
    b=[b1 b2 b3 b4 b3 b2 b1  0];

    else
        
   % 6-stage non-symmetric 5th-order: McLachla-Atela (MA92)
    a1 = 0.339839625839110000;
    a2 =-0.088601336903027329 ;
    a3 = 0.5858564768259621188 ;
    a4 =-0.603039356536491888 ;
    a5 = 0.3235807965546976394 ;
    a6 = 1-(a1+a2+a3+a4+a5); %= 0.4423637942197494587 ;

    b1 = 0.1193900292875672758 ;
    b2 = 0.6989273703824752308 ;
    b3 =-0.171312358271600775;
    b4 = 0.4012695022513534480 ;
    b5 = 0.0107050818482359840 ;
    b6 = 1-(b1+b2+b3+b4+b5); %= -0.0589796254980311632;

    a=[0  a1 a2 a3 a4 a5 a6];
    b=[b1 b2 b3 b4 b5 b6  0];

    end
elseif (s==7)  % 7-stage symmetric 6th-order ABA: Okunbor&Skeel (OS94)

	a1=-1.01308797891717472981;
	b1=0.00016600692650009894;
	a2=1.18742957373254270702;
	b2=-0.37962421426377360608;
	a3=-0.01833585209646059034;
	b3=0.68913741185181063674;
	a4=1/2-(a1+a2+a3); %=0.34399425728109261313;
	b4=1-2*(b1+b2+b3); %=0.38064159097092574080;
    a=[a1 a2 a3 a4 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b3 b2 b1  0];

elseif (s==11)   % 11-stage symmetric 6th-order BAB: Blanes&Moan (BM02)

	b1=4.146499851826237e-2;
	b2=0.198128671918067;
	b3=-4.000619210415333e-2;
	b4=7.525398430158073e-2;
	b5= -1.151138742068791e-2;
%	b6=0.236669924786931;
	b6=1/2-(b1+b2+b3+b4+b5);

	a1=0.123229775946271;
	a2= 0.290553797799558;
	a3=-0.127049212625417;
	a4=-0.246331761062075;
	a5=0.357208872795928;
%	a6=0.204777054291470;
	a6=1-2*(a1+a2+a3+a4+a5);

    a=[0  a1 a2 a3 a4 a5 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b6 b5 b4 b3 b2 b1];

elseif (s==14)   % 14-stage symmetric 6th-order ABA: Blanes&Moan (BM02)

	a1= 3.785931984061157e-2;
	a2= 0.102635633102435;
	a3=-2.586788826655869e-2;
	a4= 0.314241403071447;
	a5=-0.130144459517415;
	a6= 0.106417700369543;
	a7=-8.794243128510581e-3;
% 	a8=0.207305069056895;
	a8=1-2*(a1+a2+a3+a4+a5+a6+a7);

	b1=9.171915262446165e-2;
	b2=0.183983170005006;
	b3=-5.653436583288827e-2;
	b4=4.914688774712854e-3;
	b5=0.143761127168358;
	b6=0.328567693746804;
% 	b7=-0.196411466486454 ;
	b7=1/2-(b1+b2+b3+b4+b5+b6);

    a=[a1 a2 a3 a4 a5 a6 a7 a8 a7 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b7 b7 b6 b5 b4 b3 b2 b1  0];

elseif (s==17) 

if sa==1 
   % 17-stage symmetric 8th-order ABA: Blanes-Casas-Escorihuela (BCE24)
    a1 = 0.0520924343840339006426037968353;
    b1 = 0.145850304812644731608096609877;
    a2 = 0.225287493267702165807274831864 ;
    b2 = 0.255156544139293944162028807345;
    a3 = 0.416276189612257117795363856737 ;
    b3 = 0.0181334688208317251361460684041;
    a4 =-0.384567270213950399652168569029 ;
    b4 =-0.179040110299264554587007062749;
    a5 = 0.0997271783470514816674547589369 ;
    b5 =-0.118470801433302245053382954342;
    a6 =-0.108833834399100218757003157958 ;
    b6 = 0.186461689273821083344937258279;
    a7 = 0.222010736648991680848341975522 ;
    b7 = 0.459041581767136840219244627361;
    a8 = 0.523879522036734296002247438223 ;
    b8 =-0.003660836270318358975321459399  ;  
    
	a9=1/2-(a1+a2+a3+a4+a5+a6+a7+a8);
	b9=1-2*(b1+b2+b3+b4+b5+b6+b7+b8);

    a=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a9 a8 a7 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b7 b8 b9 b8 b7 b6 b5 b4 b3 b2 b1  0];
else
   % 17-stage symmetric 8th-order BAB: Blanes-Casas-Escorihuela (BCE24)
    
    a1 = 0.160227696073839513690970240076 ;
    b1 = 0.0514196142537210073343152693459;
    a2 = 0.306354507436867319879440957100 ;
    b2 = 0.250497030318342871458417941091;
    a3 = 0.308395508895171191756544975556 ;
    b3 = 0.512412268300327350035492806653;
    a4 = 0.120362086566233408450063177659 ;
    b4 =-0.231597138650894401279645184364;
    a5 =-0.622888687549183872072186218718 ;
    b5 = 0.116091323536875759881216298975;
    a6 = 0.635560951632990078378672016548 ;
    b6 =-0.0098365173246965763985763034283;
    a7 =-0.144226974795419229640437363913 ;
    b7 =-0.108032771466281638634277563747;
    a8 =-0.284867527074173816678992817545 ;
    b8 = 0.249039864198023642002940910070  ;  
	a9=1-2*(a1+a2+a3+a4+a5+a6+a7+a8);
	b9=1/2-(b1+b2+b3+b4+b5+b6+b7+b8);

    a=[0 a1 a2 a3 a4 a5 a6 a7 a8 a9 a8 a7 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b7 b8 b9 b9 b8 b7 b6 b5 b4 b3 b2 b1];
end

elseif (s==18)  

if sa==1 
   % 18-stage symmetric 8th-order ABA: Blanes-Casas-Escorihuela (BCE24)
    a1 = 0.0866003822712445920135805954462 ;
    b1 =-0.08;
    a2 =-0.0231572735424388070228714693753 ;
    b2 = 0.209460550048243262121199483001;
    a3 = 0.191410576083774088999564416369 ;
    b3 = 0.274887805875735483503233064415;
    a4 = 0.378895558692931579545387584925 ;
    b4 =-0.224214208870409561366168655624;
    a5 =-0.0467359566364556111599485526051; 
    b5 = 0.347657740563761656321390026010;
    a6 =-0.156198111997810415438979605642 ;
    b6 =-0.168783183866211679175007668385;
    a7 = 0.156025836895094823718831871041 ;
    b7 = 0.144209344805460873709120777707;
    a8 = 0.252844012473796333586850465807 ;
    b8 = 0.0116851121360265483381405054244;
    a9 =-0.640644212172254239866860564270  ;  
    
	a10=1-2*(a1+a2+a3+a4+a5+a6+a7+a8+a9);
	b9=1/2-(b1+b2+b3+b4+b5+b6+b7+b8);

    a=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a9 a8 a7 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b7 b8 b9 b9  b8 b7 b6 b5 b4 b3 b2 b1  0];
else
  % 18-stage symmetric 8th-order BAB: Blanes-Casas-Escorihuela (BCE24)
    
    a1 = 0.144410089394373457971755553148 ;
    b1 = 0.045;
    a2 = 0.911935520865154315536815857376 ;
    b2 = 0.459016679491512416807266107555;
    a3 =-0.00072932909837392655161199996844 ;
    b3 =-0.0456553445594333153223655352757;
    a4 =-0.930317101800698721159455541447 ;
    b4 = 0.0457031020401841003192648096559;
    a5 = 0.253804074671714046593439154323 ;
    b5 =-0.216814341025322492810152535338;
    a6 = 0.147948981530918626913598733391 ;
    b6 = 0.163168264552484857133047358600;
    a7 =-0.448814759614614928125216243784 ;
    b7 =-0.0857080319814376219389850039430;
    a8 = 0.0824123980794580106751237195418 ;
    b8 = 0.0265745810650523466142922093591;
    b9 =-0.0365538332992893220147096150675  ; 

	a9=1/2-(a1+a2+a3+a4+a5+a6+a7+a8);
	b10=1-2*(b1+b2+b3+b4+b5+b6+b7+b8+b9);

    a=[0  a1 a2 a3 a4 a5 a6 a7 a8 a9 a9 a8 a7 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b9 b8 b7 b6 b5 b4 b3 b2 b1];
end

elseif (s==19)  

if sa==1 
  % 19-stage symmetric 8th-order ABA: Blanes-Casas-Escorihuela (BCE24)
    
    a1 = 0.0505805 ;
    b1 = 0.129478606560536730662493794395;
    a2 = 0.149999 ;
    b2 = 0.222257260092671143423043559581;
    a3 =-0.0551795510771615573511026950361 ;
    b3 =-0.0577514893325147204757023246320;
    a4 = 0.423755898835337951482264998051 ;
    b4 =-0.0578312262103924910221345032763;
    a5 =-0.213495353584659048059672194633 ;
    b5 = 0.103087297437175356747933252265;
    a6 =-0.0680769774574032619111630736274 ;
    b6 =-0.140819612554090768205554103887;
    a7 = 0.227917056974013435948887201671 ;
    b7 = 0.0234462603492826276699713718626;
    a8 =-0.235373619381058906524740047732 ;
    b8 = 0.134854517356684096617882205068;
    a9 = 0.387413869179878047816794031058 ;
    b9 = 0.0287973821073779306345172160211 ;

	a10=1/2-(a1+a2+a3+a4+a5+a6+a7+a8+a9);
	b10=1-2*(b1+b2+b3+b4+b5+b6+b7+b8+b9);

    a=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a10 a9 a8 a7 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b9  b8 b7 b6 b5 b4 b3 b2 b1  0];
else   

 %17-stage symmetric 8th-order BAB: Blanes-Casas-Escorihuela (BCE24)
    a1 = 0.337548675291317241942440116575 ;
    b1 = 0.036132460472136313416730168194;
    a2 =-0.223647977575409990331768222380 ;
    b2 = 0.012697863961074113381675193011;
    a3 = 0.168949714872223740906385138015 ;
    b3 = 0.201318391240629276109068041836;
    a4 = 0.171179938816205886154783136334 ;
    b4 = 0.135683350134504233201330671671;
    a5 =-0.349765168067292877221144631312 ;
    b5 =-0.0579071833999963041504740663015;
    a6 = 0.523808861006312397712070357524 ;
    b6 =-0.0772509501792649549463874931821;
    a7 =-0.194208871063049124066394765282 ;
    b7 =-0.00264758266409925952822161203471;
    a8 =-0.323496751337931087309823477561 ;
    b8 =-0.0329844384945603065320797537355;
    a9 = 0.322817287614899749216601693799 ;
    b9 = 0.0476781560950366927530646289755 ;   

    a10=1-2*(a1+a2+a3+a4+a5+a6+a7+a8+a9);
	b10=1/2-(b1+b2+b3+b4+b5+b6+b7+b8+b9);

    a=[0  a1 a2 a3 a4 a5 a6 a7 a8 a9  a10 a9 a8 a7 a6 a5 a4 a3 a2 a1];
    b=[b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b10 b9 b8 b7 b6 b5 b4 b3 b2 b1];
end
end

