%                   coefsRKN.m
% Coefficients of selected RKN methods
% S_h= P_h^(-1) K_h P_h 
%  = (-ap_{1}, -bp_1, ..., -ap_k, -bp_k)
%    (a_{s+1}, b_s, a_s, ..., b_1, a_1)
%    (bp_k, ap_k, ..., bp_1, ap_1)
% Collected methods with s={2 4,7,11}
% S. Blanes, F. Casas and A. Murua, Splitting methods for differential equations,
%   Acta Numerica (2024)     (Section 8, Table 8.3)

function [a b ap bp] = coefsRKNproc(s,sa)
% s = Number of stages 
if (s==2)  % 2-stage non-symmetric kernel 4th-order: (BCR99)

    if sa==1
        a1 =-1/10 ;        
        a2 = 1-(a1); 
        
        b1 = 1/2-sqrt(133/132);
        b2 = 1-(b1); 
        a=[0 a1 a2 ];
        b=[b1 b2 0];
        
        z1 = 0;
        z2 = 0.5349755290809216 ;
        z3 =-0.3086327690445878 ;   
        z4 = 0.1428375011411086;
            
        y1 = 0.1937696215758170 ;
        y2 = 0.9311511462564267 ;
        y3 =-0.1053624334726687;
        y4 = 0;
        ap=-[z1 z2 z3 z4];
        bp=-[y1 y2 y3 y4];
    else
        
    end

elseif (s==4)  % 4-stage non-symmetric kernel 5th-order: (BCR99)

    if sa==1
        a1 =-3.068877096032787 ;
        a2 = 0.069136863057925 ;
        a3 = 2.132236793077397 ;    
                
        a4 = 1-(a1+a2+a3); 
        
        b1 = -0.0188 ;
        b2 =  0.228020026949214 ;
        b3 = -0.274397602327546;
        b4 = 1-(b1+b2+b3); 
            a=[0 a1 a2 a3 a4];
            b=[b1 b2 b3 b4 0];
            
        z1 = 0;
        z2 = -0.058741792287332; 
        z3 = -0.155682248638743; 
        z4 = -1.739396727797721;    
            
        y1 =-0.831482178617918 ;
        y2 = 0.869497654626095 ;
        y3 =-0.694638307539625 ;
        y4 = 0.003491184157440;
        ap=-[z1 z2 z3 z4];
        bp=-[y1 y2 y3 y4];
    else
        
    end

elseif (s==7) % 7-stage symmetric BAB kernel 6th-order: (BCR01b)
    if sa==1
        b1=0.115899400930169;
        b2=-1.21532440212000;
        b3=1.45706208067905 ;
        a1=0.244868573793901;
        a2=-0.00214552789272415;
        a3=0.301340867944477  ;     
    
	    b4=1/2-(b1+b2+b3);
	    a4=1-2*(a1+a2+a3);
        a=[0  a1 a2 a3 a4 a3 a2 a1];
        b=[b1 b2 b3 b4 b4 b3 b2 b1];
    
        z1=-0.350316247513416;
        z2=0.0744434640156453;
        z3=-0.0369370026731913 ;
        z4=-0.0597184197245884;
        z5=0.404915108936223;
        z6=-0.180941427380936 ;
        z7=-0.0346188279494959;
        z8=-(z1+z2+z3+z4+z5+z6+z7); 
        y1=0.218575120792731;
        y2=-0.370670464937763;
        y3=0.342037685653768 ;
        y4=-0.225359207496863;
        y5=0.0878524557495559;
        y6=0.195239165175742 ;
        y7=-0.155222704734044;
        y8=-(y1+y2+y3+y4+y5+y6+y7);
        ap=[z1 z2 z3 z4 z5 z6 z7 z8];
        bp=[y1 y2 y3 y4 y5 y6 y7 y8];

    else
        
    end
elseif (s==11) % 11-stage symmetric BAB kernel 6th-order: (BCR01b)

    b1=0.03906544126305366;
    b2=0.216015988434324;
    b3=-0.126717696299036 ;
    b4=-0.04128542496526060;
    b5=0.04458478096712717 ;
    a1=0.142940453575212;
    a2=0.309791505162032;
    a3=0.301210185530089 ;
    a4=-0.005822573683400349;
    a5=-0.344741324170165;
	a6=1-2*(a1+a2+a3+a4+a5); 
	b6=1/2-(b1+b2+b3+b4+b5); 
    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];
    
    z1=-0.0295940574778285;
    z2=0.0102454583206065;
    z3=0.168519324003820 ;
    z4=-0.577391651425342;
    z5=0.0991834279391326;
    z6=0.0203810695211463 ;
    z7=-0.106234446989598;
    z8=-(z1+z2+z3+z4+z5+z6+z7);  
    y1=0.175492972679660;
    y2=-0.372698829093994;
    y3=-0.00224032125918971 ;
    y4=0.0926169248899539;
    y5=-0.201446308655374;
    y6=0.216983390044259 ;
    y7=-0.0918456713646654;
    y8=-(y1+y2+y3+y4+y5+y6+y7);
    ap=[z1 z2 z3 z4 z5 z6 z7 z8];
    bp=[y1 y2 y3 y4 y5 y6 y7 y8];

    ap=[ap -ap];
    bp=[bp -bp];

end

