clear;
close all;
N = 100; %number of spartial nodes

L = 10; %x \in [-10,10]

format long;
h = 2*pi/N;
x = h*(1:N);
x = L*(x-pi)/pi; %nodes in space

kk = [0:N/2 -N/2+1:-1];
ik = 1i*kk;
ik2 = ik.*ik;

F = dftmtx(N); %FFT matrix 
iF = conj(dftmtx(N))/N; %IFFT matrix
K2 =(pi/L)^2*iF*diag(ik2)*F;
K2 = real(K2); %second derivative
K1 = (pi/L)*iF*diag(ik)*F; %first derivative
 
zero = zeros(N); 
jeden = eye(N); 

% ****** Example *******
eps = 1;
w = 1000;
M = @(t) -(1+eps*(cos(w*t))).*x.^2;
com = @(a,b) a*b - b*a ; %commutator

%Gauss-Legendre 6-th order 
x_1 = sqrt(15)/5; x_2 = 0; x_3 = -x_1;
w_1 = 5/9; w_2 = 8/9; w_3 = w_1;

c1 =0.5 - sqrt(15)/10 ; c2 = 0.5 ;c3 = 0.5 + sqrt(15)/10;
x1 = 0.08910076599011520575; x4 = 0.25995749257350879444; y1 = -0.00097618964290807330; y4 = -0.00501240016226056089; y7 = 0.00031774532164766212;
x2 = 0.24004250742649120555; x5 = 0.24789854633561981494; y2 = 0.06618969871667327349; y5 = 0.06842138031733469147;
x3 = 0.28694996084207488677; x6 = 0.00285551027560918571; y3 = 0.03862265557473451707; y6 = 0.00304401109193214959;

R1 =  @(s) iF*diag(cos((pi/L)*s*kk))*F;
R2 =  @(s) s*iF*diag(sinc((pi/L)*s*kk/pi))*F;
R3 = @(s) 1i*(K1)*iF*diag(sin((pi/L)*s*kk))*F;

load('eps01omega500_T10.mat')
Time = [1 2 4 8 16 32 64 100 150 200 400 600 800 1000 2000 5000];
error_1_2 = [];
error_2_2 = [];
error_3_2 = [];
error_4_2 = [];
error_5_2 = [];
error_6_2 = [];
error_7_2 = [];
error_8_2 = [];
error_9_2 = [];

TicToc_1_2 = [];
TicToc_2_2 = [];
TicToc_3_2 = [];
TicToc_4_2 = [];
TicToc_5_2 = [];
TicToc_6_2 = [];
TicToc_7_2 = [];
TicToc_8_2 = [];
TicToc_9_2 = [];

for l = 1:1:16
    T = Time(l);
    h = 1/T;
    hh = sqrt(h);
    
    % int_0^h M(t+s)ds
int_1 = @(s) -(x.^2.*(h*w-eps*sin(s*w) + eps* sin((h + s) *w)))/w;
% int_0^h (s-h/2)*M(t+s)ds
int_2 = @(s) -(x.^2*eps*(-2*cos(s* w)+2* cos((h + s)*w) + h*w*(sin(s*w) + sin((h + s)*w))))/(2 *w^2);
sint_1 = @(s) sqrt(int_1(s));

tic
%Strang splitting
z11= [exp(-x.^2./2)]';
z21 = zeros(1,N)';
S = [z11 ; z21];

for k = 0:T-1
    tn = k*h;
    Z1 = [diag(exp(-0.5*int_2(tn))) zero ; zero diag(exp(0.5*int_2(tn)))]*S;
    Z2 = expm([ zero h*jeden; h*K2+diag(int_1(tn)) zero])*Z1;
    S = [diag(exp(-0.5*int_2(tn))) zero ; zero diag(exp(0.5*int_2(tn)))]*Z2;
end

error_2_2 = [error_2_2  norm(Z_ref_2 - S(1:N),2)];
TicToc_2_2 = [TicToc_2_2 toc];

%BBCK 4-th
tic

c11 =0.5 - sqrt(15)/10 ; c21 = 0.5 ;c31 = 0.5 + sqrt(15)/10;

z13= [exp(-x.^2./2)]';
z23 = zeros(1,N)';

for k=0:T-1
    tn = k*h;
    W1 = diag(M(tn+c21*h));
    W2 = (sqrt(15)/3)*(diag(M(tn+c31*h))-diag(M(tn+c11*h))); 
    W3 = (10/3)*(diag(M(tn+c11*h))-2*diag(M(tn+c21*h))+diag(M(tn+c31*h)));
    
    q0 = z13; p0 = z23;
    
    q1 = q0 + h*(1/6)*p0;

    p1 = p0 + h*(0.5*K2+(0.5*W1-(1/8)*W2+(1/24)*W3))*q1;
    
    q2 = q1 + h*(2/3)*p1+(h^3)*((1/36)*K2-(7/2160)*W3)*p1;
    
    p2 = p1 + h*(0.5*K2+(0.5*W1+(1/8)*W2+(1/24)*W3))*q2;
    
    q3 = q2 + h*(1/6)*p2;
    
    z13 = q3; z23 = p2;
end 
error_3_2 = [error_3_2  norm(Z_ref_2 - z13,2)];
TicToc_3_2 = [TicToc_3_2 toc];

%Asympt^1
tic

%Gauss-Legendre , wezly, wagi
c_1 = -sqrt(3)/3; c_2 = sqrt(3)/3;

%coefficients
a0 = x.^2; a1 = eps*0.5.*x.^2; am1 = eps*0.5.*x.^2;

p00 = @(s) exp(-0.5*x.^2)*cos(s);
q00 = @(s) -exp(-0.5*x.^2)*sin(s);

q11 = @(s) -(1/1i)*(0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);
q1m1 = @(s) (1/1i)*(0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);

p10 = zeros(1,N);
q10 = zeros(1,N);

q21 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*sin(s);
q2m1 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*sin(s);

p21 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);
p2m1 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);

f = @(s) -(eps^2)/2*x.^4.*exp(-0.5*x.^2)*cos(s);
F = @(s) [ zeros(N,1);  reshape(f(s),[N,1]) ];

Om = @(s)  [diag(cos(s*sqrt(a0))) s*diag(sinc(s*sqrt(a0)/pi)) ; diag(-sin(s*sqrt(a0)).*sqrt(a0)) diag(cos(s*sqrt(a0)))]* ...
           [jeden zero; K2*s jeden];
aprox = @(s) 0.5*h*Om(0.5*h)*(Om(-0.5*h*c_1)*F(0.5*h*c_1+s+h/2)+ Om(-0.5*h*c_2)*F(0.5*h*c_2+s+h/2));
z1 = -eps*x.^2.*exp(-0.5*x.^2);
z2 = zeros(1,N);
S = [z1'; z2'];

for k =0:T-1
    tn = k*h;
    S = Om(h)*S+aprox(tn);
end
S1 = reshape(S,[1,2*N]) ;
p20 = S1(1:N);
q20 = S1(N+1:2*N);

p31 = @(s) 1/1i*eps.*x.^2.*exp(-0.5*x.^2).*sin(s);
p3m1 = @(s) -1/1i*eps.*x.^2.*exp(-0.5*x.^2).*sin(s);

p41 = @(s) eps*(4*x.^2.*exp(-0.5*x.^2).*cos(s) - exp(-0.5*x.^2).*cos(s)+0.5*x.^2.*p20);
p4m1 = @(s) eps*(4*x.^2.*exp(-0.5*x.^2).*cos(s) - exp(-0.5*x.^2).*cos(s)+0.5*x.^2.*p20);

p42 = @(s) 1/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);
p4m2 = @(s) 1/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);

ff = @(s) -eps*x.^2.*p41(s);
FF = @(s) [ zeros(N,1);  reshape(ff(s),[N,1]) ];

Om = @(s)  [diag(cos(s*sqrt(a0))) s*diag(sinc(s*sqrt(a0)/pi)) ; diag(-sin(s*sqrt(a0)).*sqrt(a0)) diag(cos(s*sqrt(a0)))]* ...
           [jeden zero; K2*s jeden];
aprox = @(s) 0.5*h*Om(0.5*h)*(Om(-0.5*h*c_1)*FF(0.5*h*c_1+s+h/2)+ Om(-0.5*h*c_2)*FF(0.5*h*c_2+s+h/2));
z11 = -(2*eps*exp(-0.5*x.^2).*(4*x.^2-1-7/16*eps*x.^4));
z21 = zeros(1,N);
S1 = [z11'; z21'];

for k = 0:T-1
    tn = k*h;
    S1 = Om(h)*S1+aprox(tn);
end
S2 = reshape(S1,[1,2*N]) ;
p40 = S2(1:N);

p51 = @(s) eps/1i*(exp(-0.5*x.^2).*sin(s).*(12*x.^2 - 4) - x.^2.*q20);
p5m1 = @(s) -eps/1i*(exp(-0.5*x.^2).*sin(s).*(12*x.^2 - 4) - x.^2.*q20);

p52 = @(s) 3/(16*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);
p5m2 = @(s) -3/(16*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);

sol_1 =  p00(1) + (1/w^2)*(p20 + 2*p21(1)*cos(1*w))+(1/w^3)*(p3m1(1)*exp(-1i*w*1)+ p31(1)*exp(1i*w*1) )+ (1/w^4)*(p40+p4m1(1)*exp(-1i*w*1)+ p41(1)*exp(1i*w*1)+p4m2(1)*exp(-1i*2*w*1)+ p42(1)*exp(1i*2*w*1))+(1/w^5)*(p5m1(1)*exp(-1i*w*1)+ p51(1)*exp(1i*w*1)+p5m2(1)*exp(-1i*2*w*1)+ p52(1)*exp(1i*2*w*1));

error_4_2 = [error_4_2  norm(Z_ref_2 - sol_1',2)];
TicToc_4_2 = [TicToc_4_2 toc];

%Asympt^2
tic

%Gauss-Legendre , wezly, wagi
c_1 = -sqrt(3)/3; c_2 = sqrt(3)/3;

%coefficients
a0 = x.^2; a1 = eps*0.5.*x.^2; am1 = eps*0.5.*x.^2;

p00 = @(s) exp(-0.5*x.^2)*cos(s);
q00 = @(s) -exp(-0.5*x.^2)*sin(s);

q11 = @(s) -(1/1i)*(0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);
q1m1 = @(s) (1/1i)*(0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);

p10 = zeros(1,N);
q10 = zeros(1,N);

q21 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*sin(s);
q2m1 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*sin(s);

p21 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);
p2m1 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);

f = @(s) -(eps^2)/2*x.^4.*exp(-0.5*x.^2)*cos(s);
F = @(s) [ zeros(N,1);  reshape(f(s),[N,1]) ];

Om = @(s)  [diag(cos(0.5*s*sqrt(a0))) diag(sin(0.5*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(0.5*s*sqrt(a0)).*sqrt(a0)) diag(cos(0.5*s*sqrt(a0)))]*[jeden zero; K2*s jeden]*[diag(cos(0.5*s*sqrt(a0))) diag(sin(0.5*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(0.5*s*sqrt(a0)).*sqrt(a0)) diag(cos(0.5*s*sqrt(a0)))];

aprox = @(s) 0.5*h*Om(0.5*h)*(Om(-0.5*h*c_1)*F(0.5*h*c_1+s+h/2)+ Om(-0.5*h*c_2)*F(0.5*h*c_2+s+h/2));
z1 = -eps*x.^2.*exp(-0.5*x.^2);
z2 = zeros(1,N);
S = [z1'; z2'];

for k =0:T-1
    tn = k*h;
    S = Om(h)*S+aprox(tn);
end
S1 = reshape(S,[1,2*N]) ;
p20 = S1(1:N);
q20 = S1(N+1:2*N);

p31 = @(s) 1/1i*eps.*x.^2.*exp(-0.5*x.^2).*sin(s);
p3m1 = @(s) -1/1i*eps.*x.^2.*exp(-0.5*x.^2).*sin(s);

q31 = @(s) 1/1i*0.5*eps.*(exp(-0.5*x.^2).*cos(s).*(-6.*x.^2+2)-x.^2.*p20);
q3m1 = @(s) -1/1i*0.5*eps.*(exp(-0.5*x.^2).*cos(s).*(-6.*x.^2+2)-x.^2.*p20);

q32 = @(s) -1/(8*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);
q3m2 = @(s) 1/(8*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);

p41 = @(s) eps*(4*x.^2.*exp(-0.5*x.^2).*cos(s) - exp(-0.5*x.^2).*cos(s)+0.5*x.^2.*p20);
p4m1 = @(s) eps*(4*x.^2.*exp(-0.5*x.^2).*cos(s) - exp(-0.5*x.^2).*cos(s)+0.5*x.^2.*p20);

p42 = @(s) 1/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);
p4m2 = @(s) 1/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);

q41 = @(s) -eps*(exp(-0.5*x.^2).*sin(s).*(-8*x.^2 + 3) + 0.5*x.^2.*q20);
q4m1 = @(s) -eps*(exp(-0.5*x.^2).*sin(s).*(-8*x.^2 + 3) + 0.5*x.^2.*q20);

q42 = @(s) 5/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);
q4m2 = @(s) 5/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);

ff = @(s) -eps*x.^2.*p41(s);
FF = @(s) [ zeros(N,1);  reshape(ff(s),[N,1]) ];

Om = @(s)  [diag(cos(0.5*s*sqrt(a0))) diag(sin(0.5*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(0.5*s*sqrt(a0)).*sqrt(a0)) diag(cos(0.5*s*sqrt(a0)))]*[jeden zero; K2*s jeden]*[diag(cos(0.5*s*sqrt(a0))) diag(sin(0.5*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(0.5*s*sqrt(a0)).*sqrt(a0)) diag(cos(0.5*s*sqrt(a0)))];
aprox = @(s) 0.5*h*Om(0.5*h)*(Om(-0.5*h*c_1)*FF(0.5*h*c_1+s+h/2)+ Om(-0.5*h*c_2)*FF(0.5*h*c_2+s+h/2));
z11 = -(2*eps*exp(-0.5*x.^2).*(4*x.^2-1-7/16*eps*x.^4));
z21 = zeros(1,N);
S1 = [z11'; z21'];

for k = 0:T-1
    tn = k*h;
    S1 = Om(h)*S1+aprox(tn);
end
S2 = reshape(S1,[1,2*N]) ;
p40 = S2(1:N);


p51 = @(s) eps/1i*(exp(-0.5*x.^2).*sin(s).*(12*x.^2 - 4) - x.^2.*q20);
p5m1 = @(s) -eps/1i*(exp(-0.5*x.^2).*sin(s).*(12*x.^2 - 4) - x.^2.*q20);

p52 = @(s) 3/(16*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);
p5m2 = @(s) -3/(16*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);

sol_2 =  p00(1) + (1/w^2)*(p20 + 2*p21(1)*cos(1*w))+(1/w^3)*(p3m1(1)*exp(-1i*w*1)+ p31(1)*exp(1i*w*1) )+ (1/w^4)*(p40+p4m1(1)*exp(-1i*w*1)+ p41(1)*exp(1i*w*1)+p4m2(1)*exp(-1i*2*w*1)+ p42(1)*exp(1i*2*w*1))+(1/w^5)*(p5m1(1)*exp(-1i*w*1)+ p51(1)*exp(1i*w*1)+p5m2(1)*exp(-1i*2*w*1)+ p52(1)*exp(1i*2*w*1));

error_5_2 = [error_5_2  norm(Z_ref_2 - sol_2',2)];
TicToc_5_2 = [TicToc_5_2 toc];

%Asympt^3
tic

%Gauss-Legendre , wezly, wagi
c_1 = -sqrt(3)/3; c_2 = sqrt(3)/3;

%Ruth coefficients
d_1_1 = 1; d_2_1 = -2/3 ; d_3_1 = 2/3;
e_1_1 = -1/24; e_2_1 = 3/4 ; e_3_1 = 7/24;

%coefficients
a0 = x.^2; a1 = eps*0.5.*x.^2; am1 = eps*0.5.*x.^2;

p00 = @(s) exp(-0.5*x.^2)*cos(s);
q00 = @(s) -exp(-0.5*x.^2)*sin(s);

q11 = @(s) -(1/1i)*(0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);
q1m1 = @(s) (1/1i)*(0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);

p10 = zeros(1,N);
q10 = zeros(1,N);

q21 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*sin(s);
q2m1 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*sin(s);

p21 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);
p2m1 =  @(s) (0.5*eps*x.^2).*exp(-0.5*x.^2)*cos(s);

f = @(s) -(eps^2)/2*x.^4.*exp(-0.5*x.^2)*cos(s);
F = @(s) [ zeros(N,1);  reshape(f(s),[N,1]) ];

Om = @(s)  [diag(cos(d_1_1*s*sqrt(a0))) diag(sin(d_1_1*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(d_1_1*s*sqrt(a0)).*sqrt(a0)) diag(cos(d_1_1*s*sqrt(a0)))]*...
           [jeden zero; e_1_1*K2*s jeden]* ...
           [diag(cos(d_2_1*s*sqrt(a0))) diag(sin(d_2_1*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(d_2_1*s*sqrt(a0)).*sqrt(a0)) diag(cos(d_2_1*s*sqrt(a0)))]*...
           [jeden zero; e_2_1*K2*s jeden]* ...
           [diag(cos(d_3_1*s*sqrt(a0))) diag(sin(d_3_1*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(d_3_1*s*sqrt(a0)).*sqrt(a0)) diag(cos(d_3_1*s*sqrt(a0)))]*...
           [jeden zero; e_3_1*K2*s jeden];
       
aprox = @(s) 0.5*h*Om(0.5*h)*(w_1*Om(-0.5*h*x_1)*F(0.5*h*x_1+s+h/2)+ ...
        w_2*Om(-0.5*h*x_2)*F(0.5*h*x_2+s+h/2)+ w_3*Om(-0.5*h*x_3)*F(0.5*h*x_3+s+h/2));
    z1 = -eps*x.^2.*exp(-0.5*x.^2);
z2 = zeros(1,N);
S = [z1'; z2'];

for k =0:T-1
    tn = k*h;
    S = Om(h)*S+aprox(tn);
end
S1 = reshape(S,[1,2*N]) ;
p20 = S1(1:N);
q20 = S1(N+1:2*N);

p31 = @(s) 1/1i*eps.*x.^2.*exp(-0.5*x.^2).*sin(s);
p3m1 = @(s) -1/1i*eps.*x.^2.*exp(-0.5*x.^2).*sin(s);

q31 = @(s) 1/1i*0.5*eps.*(exp(-0.5*x.^2).*cos(s).*(-6.*x.^2+2)-x.^2.*p20);
q3m1 = @(s) -1/1i*0.5*eps.*(exp(-0.5*x.^2).*cos(s).*(-6.*x.^2+2)-x.^2.*p20);

q32 = @(s) -1/(8*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);
q3m2 = @(s) 1/(8*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);

p41 = @(s) eps*(4*x.^2.*exp(-0.5*x.^2).*cos(s) - exp(-0.5*x.^2).*cos(s)+0.5*x.^2.*p20);
p4m1 = @(s) eps*(4*x.^2.*exp(-0.5*x.^2).*cos(s) - exp(-0.5*x.^2).*cos(s)+0.5*x.^2.*p20);

p42 = @(s) 1/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);
p4m2 = @(s) 1/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*cos(s);

q41 = @(s) -eps*(exp(-0.5*x.^2).*sin(s).*(-8*x.^2 + 3) + 0.5*x.^2.*q20);
q4m1 = @(s) -eps*(exp(-0.5*x.^2).*sin(s).*(-8*x.^2 + 3) + 0.5*x.^2.*q20);

q42 = @(s) 5/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);
q4m2 = @(s) 5/(16)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);

ff = @(s) -eps*x.^2.*p41(s);
FF = @(s) [ zeros(N,1);  reshape(ff(s),[N,1]) ];

Om = @(s)  [diag(cos(d_1_1*s*sqrt(a0))) diag(sin(d_1_1*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(d_1_1*s*sqrt(a0)).*sqrt(a0)) diag(cos(d_1_1*s*sqrt(a0)))]*...
           [jeden zero; e_1_1*K2*s jeden]* ...
           [diag(cos(d_2_1*s*sqrt(a0))) diag(sin(d_2_1*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(d_2_1*s*sqrt(a0)).*sqrt(a0)) diag(cos(d_2_1*s*sqrt(a0)))]*...
           [jeden zero; e_2_1*K2*s jeden]* ...
           [diag(cos(d_3_1*s*sqrt(a0))) diag(sin(d_3_1*s*sqrt(a0))./sqrt(a0)) ; diag(-sin(d_3_1*s*sqrt(a0)).*sqrt(a0)) diag(cos(d_3_1*s*sqrt(a0)))]*...
           [jeden zero; e_3_1*K2*s jeden];
       
aprox = @(s) 0.5*h*Om(0.5*h)*(w_1*Om(-0.5*h*x_1)*FF(0.5*h*x_1+s+h/2)+ ...
        w_2*Om(-0.5*h*x_2)*FF(0.5*h*x_2+s+h/2)+ w_3*Om(-0.5*h*x_3)*FF(0.5*h*x_3+s+h/2));
    
    z11 = -(2*eps*exp(-0.5*x.^2).*(4*x.^2-1-7/16*eps*x.^4));
z21 = zeros(1,N);
S1 = [z11'; z21'];

for k = 0:T-1
    tn = k*h;
    S1 = Om(h)*S1+aprox(tn);
end
S2 = reshape(S1,[1,2*N]) ;
p40 = S2(1:N);

p51 = @(s) eps/1i*(exp(-0.5*x.^2).*sin(s).*(12*x.^2 - 4) - x.^2.*q20);
p5m1 = @(s) -eps/1i*(exp(-0.5*x.^2).*sin(s).*(12*x.^2 - 4) - x.^2.*q20);

p52 = @(s) 3/(16*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);
p5m2 = @(s) -3/(16*1i)*eps^2*x.^4.*exp(-0.5*x.^2).*sin(s);

sol_3 =  p00(1) + (1/w^2)*(p20 + 2*p21(1)*cos(1*w))+(1/w^3)*(p3m1(1)*exp(-1i*w*1)+ p31(1)*exp(1i*w*1) )+ (1/w^4)*(p40+p4m1(1)*exp(-1i*w*1)+ p41(1)*exp(1i*w*1)+p4m2(1)*exp(-1i*2*w*1)+ p42(1)*exp(1i*2*w*1))+(1/w^5)*(p5m1(1)*exp(-1i*w*1)+ p51(1)*exp(1i*w*1)+p5m2(1)*exp(-1i*2*w*1)+ p52(1)*exp(1i*2*w*1));

error_6_2 = [error_6_2  norm(Z_ref_2 - sol_3',2)];
TicToc_6_2 = [TicToc_6_2 toc];

% Xi^3 
    
R1h = R1(h);
R2h = R2(h);
R3h = R3(h);

z_pop= [exp(-x.^2./2)]';
z_pop_der = zeros(1,N)';
z_zapop = zeros(1,N)';
S7 = [z_pop;z_pop_der];

int_1 = @(t) -x.^2.*(h + (eps*(sin(w*(h + t)) - sin(t*w)))/w);
int_2 = @(t) -x.^2.*( h^2/2 + eps*((cos(w*(h + t)) + h*w*sin(w*(h + t)))/w^2 - cos(t*w)/w^2));
int_3 = @(t) -x.^2.*( h^3/3 + eps*((h^2*w^2*sin(w*(h + t)) - 2*sin(w*(h + t)) + 2*h*w*cos(w*(h + t)))/w^3 + (2*sin(t*w))/w^3));

tn = 0;
S1 = [ R1h*z_pop+R2h*z_pop_der + (h*ifft(diag(sinc((pi/L)*h*kk/pi))*fft(diag(int_1(tn))))-ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_2(tn))))+(1/(2*h))*(-1*diag(int_3(tn)) +ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_3(tn))))))*z_pop+...
                                     (h*ifft(diag(sinc((pi/L)*h*kk/pi))*fft(diag(int_2(tn))))+(1/(2*h))*(-h*diag(int_3(tn)) - h*ifft(diag(sinc((pi/L)*h*kk/pi))*fft(diag(int_3(tn))))))*z_pop_der+... 
                                     ((1/(2*h))*(-1*0.5*h^2)*diag(int_3(tn)))*(K2*z_pop);...
           R3h*z_pop+R1h*z_pop_der + (ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_1(tn))))-1i*K1*ifft(diag(sin((pi/L)*h*kk))*fft(diag(int_2(tn))))+(1/(2*h))*(1i*K1*ifft(diag(sin((pi/L)*h*kk))*fft(diag(int_3(tn))))))*z_pop+...
                                     (ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_2(tn))))+(1/(2*h))*(1*diag(int_3(tn)) - ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_3(tn))))))*z_pop_der+... 
                                     ((1/(2*h))*h*ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_3(tn)))))*(K2*z_pop)];
z_zapop = z_pop_der;
z_pop = S1(1:N); z_pop_der = S1(N+1:2*N);
S7 = [z_pop;z_pop_der];

tic
for k=1:T-1
    tn = k*h;
    S7 = [ R1h*z_pop+R2h*z_pop_der + (h*ifft(diag(sinc((pi/L)*h*kk/pi))*fft(diag(int_1(tn))))-ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_2(tn))))+(1/(2*h))*(-1*diag(int_3(tn)) +ifft(diag(cos(h*kk))*fft(diag(int_3(tn))))))*z_pop+...
                                     (h*ifft(diag(sinc((pi/L)*h*kk/pi))*fft(diag(int_2(tn))))+(1/(2*h))*(-h*diag(int_3(tn)) - h*ifft(diag(sinc((pi/L)*h*kk/pi))*fft(diag(int_3(tn))))))*z_pop_der+... 
                                     ((1/(2*h))*(-1*0.5*h^2)*diag(int_3(tn)))*((1/(h))*(z_pop_der - z_zapop));...
           R3h*z_pop+R1h*z_pop_der + (ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_1(tn))))-1i*K1*ifft(diag(sin((pi/L)*h*kk))*fft(diag(int_2(tn))))+(1/(2*h))*(1i*K1*ifft(diag(sin((pi/L)*h*kk))*fft(diag(int_3(tn))))))*z_pop+...
                                     (ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_2(tn))))+(1/(2*h))*(1*diag(int_3(tn)) - ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_3(tn))))))*z_pop_der+... 
                                     ((1/(2*h))*h*ifft(diag(cos((pi/L)*h*kk))*fft(diag(int_3(tn)))))*(1/(h))*(z_pop_der - z_zapop)];
    z_zapop = z_pop_der;
     z_pop = S7(1:N); z_pop_der = S7(N+1:2*N);      
end
error_7_2 = [error_7_2  norm(Z_ref_2 - z_pop,2)];
TicToc_7_2 = [TicToc_7_2 toc];

% int_0^h M(t+s)ds
int_1 = @(s) -(x.^2.*(h*w-eps*sin(s*w) + eps* sin((h + s) *w)))/w;
% int_0^h (s-h/2)*M(t+s)ds
int_2 = @(s) -(x.^2*eps*(-2*cos(s* w)+2* cos((h + s)*w) + h*w*(sin(s*w) + sin((h + s)*w))))/(2 *w^2);
sint_1 = @(s) sqrt(int_1(s));

%Gamma_2
tic

z1= [exp(-x.^2./2)]';
z2 = zeros(1,N)';

for k= 0:T-1
    tn = k*h; 
    q0 = diag(exp(-0.5*int_2(tn)))*z1; p0 = diag(exp(0.5*int_2(tn)))*z2;
    
    p1 = (1/6)*h*K2*q0+p0;
    
    q1 = diag(cosh(hh*sint_1(tn)./2))*q0 + diag(hh*sinh(hh*sint_1(tn)/2)./(sint_1(tn)))*p1;
    p2 = diag(cosh(hh*sint_1(tn)./2))*p1 + diag(sint_1(tn).*sinh(hh*sint_1(tn)/2)./hh)*q0;
    
    p3 = ((2/3)*h*K2+(2/72)*h*h*h*K2*K2)*q1+p2;
    
    q2 = diag(cosh(hh*sint_1(tn)./2))*q1 + diag(hh*sinh(hh*sint_1(tn)/2)./(sint_1(tn)))*p3;
    p4 = diag(cosh(hh*sint_1(tn)./2))*p3 + diag(sint_1(tn).*sinh(hh*sint_1(tn)/2)./hh)*q1;
    
    p5 = (1/6)*h*K2*q2+p4;
    
    q3 = diag(exp(-0.5*int_2(tn)))*q2; p6 = diag(exp(0.5*int_2(tn)))*p5;
    
    z1 = q3;
    z2 = p6;
end

error_8_2 = [error_8_2  norm(Z_ref_2 - z1,2)];
TicToc_8_2 = [TicToc_8_2 toc];

tic
%Gamma_1
z11= [exp(-x.^2./2)]';
z21 = zeros(1,N)';

for k= 0:T-1
    tn = k*h;  
    q0 = diag(exp(-0.5*int_2(tn)))*z11; p0 = diag(exp(0.5*int_2(tn)))*z21;
    
    p1 = 1/6*(h*K2+diag(int_1(tn)))*q0 + p0;
    q1 = q0+0.5*h*jeden*p1;
    
    p2 = ((2/3)*(h*K2+diag(int_1(tn)))+(2/72)*h*(h*K2+diag(int_1(tn)))^2)*q1+p1;
    q2 = q1+0.5*h*jeden*p2;
    
    p3 = 1/6*(h*K2+diag(int_1(tn)))*q2 + p2;
    
    q3 = diag(exp(-0.5*int_2(tn)))*q2; 
    
    p4 = diag(exp(0.5*int_2(tn)))*p3;
    
    z11 = q3;
    z21 = p4;
end

error_9_2 = [error_9_2  norm(Z_ref_2 - z11,2)];
TicToc_9_2 = [TicToc_9_2 toc];

end

% figure
% subplot(1,2,1)
% loglog(Time.^(-1), error_2_2,'-o',Time.^(-1), error_3_2,'-o',Time.^(-1), error_4_2,'-v',Time.^(-1), error_5_2,'-v',Time.^(-1), error_6_2,'-v',Time.^(-1), error_7_2,'-*',Time.^(-1), error_8_2,'-*',Time.^(-1), error_9_2,'-p',Time.^(-1),Time.^(-2),'k-.',Time.^(-1),Time.^(-3),'k--',Time.^(-1),Time.^(-4),'k:','markersize',6), grid on, xlabel ' Time step h', ylabel 'error', title ' \epsilon = 1 , \omega = 1000 '
% legend('Strang','BBCK^{[4]}','Asymp^{[1]}','Asymp^{[2]}','Asymp^{[3]}', '\Xi^{[3]}','\Gamma^{[4]}_2','\Gamma^{[4]}_1','O(h^2)','O(h^3)','O(h^4)')
% subplot(1,2,2)
% loglog(TicToc_2_2, error_2_2,'-o',TicToc_3_2, error_3_2,'-o',TicToc_4_2, error_4_2,'-v',TicToc_5_2, error_5_2,'-v',TicToc_6_2, error_6_2,'-v',TicToc_7_2, error_7_2,'-*',TicToc_8_2, error_8_2,'-*',TicToc_9_2, error_9_2,'-p','markersize',6), grid on, xlabel ' Seconds', ylabel 'error', title ' \epsilon = 1 , \omega = 1000 '

save('wave_200_error_2_2','error_2_2');
save('wave_200_error_3_2','error_3_2');
save('wave_200_error_4_2','error_4_2');
save('wave_200_error_5_2','error_5_2');
save('wave_200_error_6_2','error_6_2');
save('wave_200_error_7_2','error_7_2');
save('wave_200_error_8_2','error_8_2');
save('wave_200_error_9_2','error_9_2');

save('wave_200_tictoc_2_2','TicToc_2_2');
save('wave_200_tictoc_3_2','TicToc_3_2');
save('wave_200_tictoc_4_2','TicToc_4_2');
save('wave_200_tictoc_5_2','TicToc_5_2');
save('wave_200_tictoc_6_2','TicToc_6_2');
save('wave_200_tictoc_7_2','TicToc_7_2');
save('wave_200_tictoc_8_2','TicToc_8_2');
save('wave_200_tictoc_9_2','TicToc_9_2');