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;
M2 = @(t) -(1+eps*(cos(w*t))).*x.^2;
com = @(a,b) a*b - b*a ; %commutator

% %Gauss-Legendre , wezly, wagi
% x_1 = sqrt(15)/5; x_2 = 0; x_3 = -x_1;
% w_1 = 5/9; w_2 = 8/9; w_3 = w_1;

%reference solution 
TR = 1000000;
h = 1/TR;

A = @(s) [ zero jeden; K2+diag(M2(s)) zero];

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


for k=0:TR-1
    tn = k*h;
    A1 = A(tn+(0.5-sqrt(15)/10)*h); A2 = A(tn+0.5*h); A3 = A(tn+(0.5+sqrt(15)/10)*h);
    alfa1 = h*A2;
    alfa2 = sqrt(15)/3*h*(A3-A1);
    alfa3 = 10*h/3*(A3-2*A2+A1);
    C1 = com(alfa1,alfa2);
    C2 = -1/60*com(alfa1,2*alfa3+C1);
    S = expm(alfa1+1/12*alfa3+(1/240)*com(-20*alfa1-alfa3+C1,alfa2+C2))*S;
end
Z_ref_2 = S(1:N);
save('eps01omega500_T10','Z_ref_2')