%                           coefsNI.m
% Coefficients of selected splitting methods for near-integrable systems
%   y'=f_0(y) + eps f_I(y) = A + eps B
%  S_h = psi^B_{b_{s+1} h} psi^A_{a_{s+1} h}···psi^B_{b_1 h} psi^A_{a_1 h}
% Methods of order (s1,s2,s3) ~ eps h^(s1) + eps^2 h^(s2) + eps^3 h^(s3)
% The input values s1, s2 suffice to label the methods
% ab=1: ABA methods
% ab=0: BAB methods
% (s1,s2)={(4,2),(6,2),(8,2),(10,2),(6,4),(8,4),(10,4),(8,6),(10,6)} 
% Output: the coefficients a_i, b_i, and th enumber of stages, s
% S. Blanes, F. Casas and A. Murua, Splitting methods for differential equations,
%   Acta Numerica (2024)     (Section 8, Table 8.5)

 function [a b s] = coefsNI(s1,s2,ab)
% s = Number of stages 
if (s1==2 & s2==2)  
        s=1;
    if (ab==1) % Strang ABA (WH91)
       a=[1/2 , 1/2];
       b=[1 , 0];
    else  % Strang BAB
       a=[0, 1 ];
       b=[1/2 , 1/2];
    end
elseif (s1==4 & s2==2) 
        s=2;
    if (ab==1) % 2-stage (4,2) ABA method (MCL95a)
        a1=(3-sqrt(3))/6;
        a2=1-2*a1;
        a3=a1;
        b1=1/2;
        b2=b1;
       a=[a1 , a2 , a3];
       b=[b1 , b2 , 0];
    else  % 2-stage (4,2) BAB method (MCL95a)
        a1=1/2;
        a2=a1;
        b1=1/6;
        b2=1-2*b1;
        b3=b1;
       a=[0, a1 , a2 ];
       b=[b1 , b2 , b3];
    end
elseif (s1==6 & s2==2)  
        s=3;
    if (ab==1) % 3-stage (6,2) ABA method (MCL95a)
        a1=1/2-sqrt(3/20);
        a2=1/2-a1;
        a3=a2;
        a4=a1;
        b1=5/18;
        b2=1-2*b1;
        b3=b1;
       a=[a1 , a2 , a3 , a4];
       b=[b1 , b2 , b3 , 0];
    else  % 3-stage (6,2) BAB method (MCL95a)
        a1=(5-sqrt(5))/10;
        a2=1-2*a1;
        a3=a1;
        b1=1/12;
        b2=1/2-b1;
        b3=b2;
        b4=b1;
       a=[0, a1 , a2 , a3 ];
       b=[b1 , b2 , b3 , b4];
    end
elseif (s1==8 & s2==2)  
        s=4;
    if (ab==1) % 4-stage (8,2) ABA method (MCL95a)
          w1=sqrt((15+2*sqrt(30))/35);   w2=sqrt((15-2*sqrt(30))/35); 
        a1=(1-w1)/2;
        a2=(w1-w2)/2;
        a3=1-2*(a1+a2);
        a4=a2;
        a5=a1;
        b1=(18-sqrt(30))/72;
        b2=1/2-b1;
        b3=b2;
        b4=b1;
       a=[a1 , a2 , a3 , a4 , a5];
       b=[b1 , b2 , b3 , b4 , 0];
    else % 4-stage (8,2) BAB method (MCL95a)
        a1=1/2 - sqrt(3/28);
        a2=1/2-a1;
        a3=a2;
        a4=a1;
        b1=1/20;
        b2=49/180;
        b3=1-2*(b1+b2);
        b4=b2;
        b5=b1;
       a=[0, a1 , a2 , a3 , a4 ];
       b=[b1 , b2 , b3 , b4 , b5];
    end
elseif (s1==8 & s2==4)  
        s=5;
    if (ab==1) % 5-stage (8,4) ABA method (MCL95a)
        a1= 0.07534696026989288842;
        a2= 0.51791685468825678230;
        a3=1/2-(a1+a2);
        a4=a3;
        a5=a2;
        a6=a1;
        b1= 0.19022593937367661925;
        b2= 0.84652407044352625706;
        b3=1-2*(b1+b2);
        b4=b2;
        b5=b1;
       a=[a1 , a2 , a3 , a4 , a5 , a6];
       b=[b1 , b2 , b3 , b4 , b5 , 0];
    else % 5-stage (8,4) BAB method (MCL95a)
        a1=-0.00758691311877447385;
        a2=0.31721827797316981388;
        a3=1-2*(a1+a2);
        a4=a2;
        a5=a1;
        b1=0.81186273854451628884;
        b2=-0.67748039953216912289;
        b3=1/2-(b1+b2);
        b4=b3;
        b5=b2;
        b6=b1;
       a=[0, a1 , a2 , a3 , a4 , a5 ];
       b=[b1 , b2 , b3 , b4 , b5 , b6];
    end
elseif (s1==10 & s2==4)   % 7-stage (10,4) ABA method (BCF13)
	imet=7; iAB=1; s=7;
        a1=0.04706710064597250612947887637243678556564;
        b1= 0.1188819173681970199453503950853885936957;
        a2= 0.1847569354170881069247376193702560968574;
        b2=0.2410504605515015657441667865901651105675 ;
        a3= 0.2827060056798362053243616565541452479160 ;
        b3= -0.2732866667053238060543113981664559460630 ;
       a4=1/2-(a1+a2+a3) ;
%         a4= -0.01453004174289681837857815229683813033908;
       b4=  1-2*(b1+b2+b3);
%         b4=0.8267085775712504407295884329818044835997;
       a=[a1 , a2 , a3 , a4 , a4 , a3 , a2 , a1];
       b=[b1 , b2 , b3 , b4 , b3 , b2 , b1 ,  0];

elseif (s1==8 & s2==6)  
        s=7;
    if (ab==1)  % 7-stage (8,6,4) ABA method (BCF13)
        a1 = 0.0711334264982231177779387300061549964174;
        a2 = 0.241153427956640098736487795326289649618;
        a3 = 0.521411761772814789212136078067994229991;
        % a4 =-0.333698616227678005726562603400438876027;
        b1 = 0.183083687472197221961703757166430291072;
        b2 = 0.310782859898574869507522291054262796375;
        b3 =-0.0265646185119588006972121379164987592663;
        % b4 = 0.0653961422823734184559721793911134363710  ;      
        a4=1/2-(a1+a2+a3);
        b4=1-2*(b1+b2+b3);
       a=[a1 , a2 , a3 , a4 , a4 ,a3 , a2 , a1];
       b=[b1 , b2 , b3 , b4 , b3 , b2 , b1 , 0];
    else

    end
elseif (s1==10 & s2==6)  
        s=8;
    if (ab==1)  % 8-stage (10,6,4) ABA method (BCF13)
        a1= 0.03809449742241219545697532230863756534060;
        a2= 0.1452987161169137492940200726606637497442;
        a3= 0.2076276957255412507162056113249882065158;
        a4= 0.435909703651526159223154862401065184400;
        a5=1-2*(a1+a2+a3+a4);
        b1= 0.09585888083707521061077150377145884776921;
        b2= 0.2044461531429987806805077839164344779763;
        b3= 0.2170703479789911017143385924306336714532;
        b4=1/2-(b1+b2+b3);
       a=[a1 , a2 , a3 , a4 , a5 , a4 , a3 , a2 , a1];
       b=[b1 , b2 , b3 , b4 , b4 ,  b3 , b2 , b1 , 0];
    else

    end
end

