% Berechnet die Losung der 1D Schroedingergleichung
%
% d_t v(t,x) = i d_xx v(t,x) + i*Pot(x) * v(t,x) 
%     v(0,x) = v_0(x)

T = 0.2;   % Endzeitpunkt
N = 512; dx = 2*pi/N; x = dx*(1:N)'; 

% Startwerte
v_0  = max(0,1-abs(x-pi)/2);  % Hutfunktion
%v_0  = exp(-3*(x-pi).^2);     % Gaussverteilung

column = [0 .5*(-1).^(1:N-1).*cot((1:N-1)*dx/2)]';
D = toeplitz(column,column([1 N:-1:2])); % 1. Ableitung
D2 = D*D; % 2. Ableitung

PotV = sin(x);
PotM = diag(PotV);

% Referenzloesung v(T) = exp( T *(i d_xx + i M_pot) )v_0
vref = expm(1i*(D2+PotM)*T)*v_0;

errorT = []; % Fehler fuer Trotterprodukt
errorS = []; % Fehler fuer Strangsplitting

dtvek = [];

for Nt = [20 50 80 100 150 200 300 500 800 1000 1200 1600 2000]
	
	dt = T/Nt;
	vT = v_0;
	vS = v_0;
	
	for tstep = 1:Nt
		%% Trotter Splitting
		%T(h)
		w_hat = exp(1i*dt*(1i*[0:N/2-1 0 -N/2+1:-1]').^2).*fft(vT);
		vT = ifft(w_hat);
		%S(h)
		vT = exp(1i*dt*PotV).*vT;
		
		%% Strang splitting
		%T(h/2)
		v_hat = exp(1i*dt/2*(1i*[0:N/2-1 0 -N/2+1:-1]').^2).*fft(vS);
		vS = ifft(v_hat);
		%S(h)
		vS = exp(1i*dt*PotV).*vS;
		%T(h/2)
		w_hat = exp(1i*dt/2*(1i*[0:N/2-1 0 -N/2+1:-1]').^2).*fft(vS);
		vS = ifft(w_hat);
		
	end
	errorT = [errorT sqrt((vT-vref)'*(vT-vref))];
	errorS = [errorS sqrt((vS-vref)'*(vS-vref))];
	dtvek = [dtvek dt];
end

figure(1)
loglog(dtvek,errorT,'bx-',dtvek,errorS,'rx-');
xlabel('dt'); ylabel('error'); legend('Trotter','Strang');
shg