N = 256;
T = 1;
l = -N : N-1;
x = l*pi/N;
V = 1-cos(x);

% arrange initial values for the ODE
% ----------------------------------
r = rand(1,2*N);
v = r/norm( r,2);
U01 = ( diag(1./sqrt(l.^2+1) ) * v' ) / norm( diag(1./sqrt(l.^2+1) ) * v' , 2);%red%
U02 = ( diag(1./(l.^2+1) ) * v' ) / norm( diag(1./(l.^2+1) ) * v' , 2);%blue% 

% exact solution of the ODE with initial values U01 and U02
% ---------------------------------------------------------
U_exact1 = @(t) expm( ( diag(-l.^2) - fft(eye(2*N)) * diag(V) * ifft(eye(2*N)) )*t*i )*U01;
U_exact2 = @(t) expm( ( diag(-l.^2) - fft(eye(2*N)) * diag(V) * ifft(eye(2*N)) )*t*i )*U02;

% Stepsizes
% ---------
stepsize = [1,1/3, 1/3^2, 1/3^3, 1/3^4, 1/3^5, 1/3^6];

% U(:,k+1) is solution in t=k*tau
% abserr(k) is error for stepsize(k) 

for k = 1 : length(stepsize)
% compute error in t=1 (Marchuk-Strang splitting) 
% -----------------------------------------------
[u1_MS,U1_MS] = solver(T,N,V,U01,stepsize(k),2); 
abserr1_MS(k) = max(max(abs( U_exact1(stepsize(k)) - U1_MS(:,2) )));
[u2_MS,U2_MS] = solver(T,N,V,U02,stepsize(k),2); 
abserr2_MS(k) = max(max(abs( U_exact2(stepsize(k)) - U2_MS(:,2) )));
end
% plot error in t=1 (Marchuk-Strang splitting)
% --------------------------------------------
figure;
hold on
title('Marchuk-Strang splitting')
plot(stepsize,abserr1_MS,'r-+')
plot(stepsize,abserr2_MS,'b-o')
xlabel('tau');
ylabel('error');
set(gca, 'color', 'white')
set(gca, 'XScale', 'log')
set(gca, 'yScale', 'log')

plot(stepsize,abserr1_MS./stepsize,'r--+')
plot(stepsize,abserr2_MS./stepsize,'b--o')