%                   coefsS.m
% Coefficients of selected Method-Adjoint compositions
%  S_h = chi_{ga_{2s}h}·chi^*_{ga_{2s-1}h}··· chi_{ga_{2}h}·chi^*_{ga_{1}h}
% Methods with number of stages s={1 4,5,6,10}
% S. Blanes, F. Casas and A. Murua, Splitting methods for differential equations,
%   Acta Numerica (2024)     (Section 8, Table 8.2)

function [ga] = coefsS(s)
% s = Number of stages 
if (s==1)  % 2nd-order: Strang
   ga=[1/2 1/2];
elseif (s==2)  % 2nd-order: McLachlan (MCL95b/OMF03)
    y=(2*sqrt(326)-36)^(1/3); 
    b1=(y^2+6*y-2)/(12*y); 
  b=[b1,1-2*(b1),b1];
    a1=1/2; 
  a=[a1,a1];
    ga(1)=b(1);
    ga(2)=a(1)-ga(1);
    for i=2:2;
        ga(2*i-1)=b(i)-ga(2*i-2);
        ga(2*i)=a(i)-ga(2*i-1);
    end

elseif (s==4)  % 4th-order: McLachlan (MCL95b)
    b1=(642+sqrt(471))/3924; 
    b2=121*(12-sqrt(471))/3924;
  b=[b1,b2,1-2*(b1+b2),b2,b1];
    a1=6/11; 
  a=[a1,1/2-a1,1/2-a1,a1];
    ga(1)=b(1);
    ga(2)=a(1)-ga(1);
    for i=2:4;
        ga(2*i-1)=b(i)-ga(2*i-2);
        ga(2*i)=a(i)-ga(2*i-1);
    end

elseif (s==5)   % 4th-order: McLachlan (MCL95b)
  nstages=5;
    b1=(14-sqrt(19))/108; 
    b2=(20-7*sqrt(19))/108;
  b=[b1,b2,1/2-(b1+b2),1/2-(b1+b2),b2,b1];
    a1=2/5; 
    a2=-1/10; 
  a=[a1,a2,1-2*(a1+a2),a2,a1];
    ga(1)=b(1);
    ga(2)=a(1)-ga(1);
    for i=2:5;
        ga(2*i-1)=b(i)-ga(2*i-2);
        ga(2*i)=a(i)-ga(2*i-1);
    end

elseif (s==6)   % 4th-order: Blanes&Moan (BM02)
    a1=0.209515106613362; 
    a2=-0.143851773179818;
  a=[a1,a2,1/2-(a1+a2),1/2-(a1+a2),a2,a1];
    b1=0.07920369643119565; 
    b2=0.353172906049774; 
    b3=-0.04206508035771952;
  b=[b1,b2,b3,1-2*(b1+b2+b3),b3,b2,b1];
    ga(1)=b(1);
    ga(2)=a(1)-ga(1);
    for i=2:6;
        ga(2*i-1)=b(i)-ga(2*i-2);
        ga(2*i)=a(i)-ga(2*i-1);
    end

elseif (s==10)   % 6th-order: Blanes&Moan (BM02)
    a1=0.148816447901042; 
    a2=-0.132385865767784;
    a3=6.730760469218501e-2; 
    a4=0.432666402578175;
%     a5=-1.640458940361793d-002;
  a=[a1,a2,a3,a4,1/2-(a1+a2+a3+a4),1/2-(a1+a2+a3+a4),a4,a3,a2,a1];
    b1=5.026276440039223e-2; 
    b2=0.413514300428344; 
    b3=4.507988979439766e-2;
    b4=-0.188054853819569;
    b5=0.541960678450780;
%     b6=-0.72552555850868;
  b=[b1,b2,b3,b4,b5,1-2*(b1+b2+b3+b4+b5),b5,b4,b3,b2,b1];
    ga(1)=b(1);
    ga(2)=a(1)-ga(1);
    for i=2:10;
        ga(2*i-1)=b(i)-ga(2*i-2);
        ga(2*i)=a(i)-ga(2*i-1);
    end
end

