function [ solutions, error, relerror, relerrorsquare ] = schroedinger_operatorsplitting(N,h,t,scheme,u0)
%SCHROEDINGER_OPERATORSPLITTING Operator splitting for Schroedinger
%equation

% Solves the space-discretizated Schroedinger equation approximately via 
% operator splitting.

% N = Number of mesh points,
% h = time steps,
% t = time point of interest,
% scheme = operator splitting scheme, 1 for sequation, 2 for Marchuk-Strang
% u0 = initial value, column vector of length 2N

% Set time interval
a = pi;
x = zeros(2*N,1);
for k=1:2*N
    x(k) = (-N-1+k)*a/N;
end

% Calculate matrices
D = i*diag([-N:(N-1)]);
Dsquare = D*D;
% Potential V(x) = 1-cos(x)
%V = diag(x+*pi);
V = diag(1-cos(x));


% Splitting
n = floor(t/h);
t = h*n;

% Give ALL the solutions
solutions = zeros(2*N,n+1);
solutions(:,1) = u0;

uold = u0;
% Sequential splitting
if (scheme == 1) 
    % Matrix exponentials
    eD2 = expm(i*h*Dsquare);
    eV = expm(-i*h*V);
    for k=1:n
        unew = fft(eV*ifft(eD2*uold));
        uold = unew;
    end
% Marchuk-Strang splitting
else
    % Matrix exponentials
    eD2h = expm(i*(h/2)*Dsquare);
    eV = expm(-i*h*V);
    for k=1:n
        solutions(:,k+1) = eD2h*fft(eV*ifft(eD2h*solutions(:,k)));
    end
end

if (nargout > 1)
    % Calculate correct solution of space-discretizated ODE
    uright = expm(i*t*(Dsquare - dftmtx(2*N)*V*conj(dftmtx(2*N))/(2*N)))*u0;

    % Error estimate
    error = norm(solutions(:,n+1)-uright);
    relerror = error/h;
    relerrorsquare = relerror/h;
end
