function aufgabe9a(N)

if (nargin < 1)
    N = 256;
end

V = @(x) -cos(x/2);
%s = 0.1;
psi0 = @(x) exp(-10*(x-0.5).^2) * (20/pi)^(1/4);
mu = 1/0.1^2;

k = (-N/2):(N/2-1);
x = k * 2 * pi / N;
w = ones(size(x)) * (x(2) - x(1));
w(1) = w(1)/2;
w(end) = w(end)/2;

plot(x, V(x), '-k', x, psi0(x), '-r'); xlim([-pi pi]); drawnow;
fprintf('norm is: %g\n', w * psi0(x').^2);

DK = 1/(2 * mu) * k.^2;

tic
F = 1/N * exp(-2i*pi/N * k' * k);
FI = exp(2i*pi/N * k' * k);
toc

%tic;
%V1 = F * diag(V(x)) * FI;
%toc

x1 = linspace(-pi, pi, N+1);
w1 = ones(size(x1)) * (x1(2) - x1(1));
w1(1) = w1(1)/2;
w1(end) = w1(end)/2;
tic
V2 = 1/(2*pi) * ones(N,1) * (w1 .* V(x1)) .* ...
    exp(-1i*k'*x1) * transpose(exp(1i*k'*x1));
toc

c0 = F * psi0(x');

%c1 = 1/(2*pi) * exp() (w1 * psi0(x))'

%f = @(t, y) -1i * (diag(DK) + V1) * y;
%tic
%[t1, ct1] = ode45(f, [0 100], c0);
%toc
%tic
%[t2, ct2] = ode15s(f, [0 100], c0);
%toc

k = [0:(N/2-1) (-N/2):-1];
DK = 1/(2 * mu) * k.^2;
c0 = fft(psi0(x'));
f = @(t, y) -1i * (DK' .* y + fft(V(x1(2:end))' .* ifft(y)));

fprintf('V=cos(x), ode45 ');
tic
[t3, ct3] = ode45(f, [0 100], c0);
toc
view2D(1, ct3, t3, x);
title('ode45');

fprintf('V=cos(x), ode15s ');
tic
[t4, ct4] = ode15s(f, [0 100], c0);
toc
view2D(2, ct4, t4, x);
title('ode15s');

h = @(x, a) 0.5*(sign(x-a)+1);
hm1 = @(x, a) -0.5*(sign(x-a)-1);
%Vs = @(x) hm1(x,-5).*((x+5).^2 + 2) + h(x,-5).*h(x,0).*ceil(-2/5*x) ...
%     + (max(x,5)-5).^2;

%V = @(x) Vs(x*pi/10);
V = @(x) abs(x);
%plot(x, V(x));

c0 = fft(psi0(x'));
f = @(t, y) -1i * (DK' .* y + fft(V(x1(2:end))' .* ifft(y)));
fprintf('V=|x|, ode45s ');
tic
[t5, ct5] = ode45(f, [0 100], c0);
toc
view2D(3, ct5, t5, x);
title('V(x) = |x|');

figure(4); clf;
Z = ifft(transpose(ct5));
NX = length(x); IX = 1:ceil(NX/20):NX;
NT = length(t5); IT = 1:ceil(NT/50):NT;
[X T] = ndgrid(x(IX), t5(IT));
quiver(X, T, real(Z(IX,IT)), imag(Z(IX,IT)));

return

function view2D(fig, ct, t, x)
[X T] = ndgrid(x, t);
Z = ifft(transpose(ct));
C = thaller(Z, 0.5);

hf = figure(fig); clf;
contour3(X, T, abs(Z), 4, '-w'); % Konturlinien
hs = surface(X, T, abs(Z), hsv2rgb(C), 'EdgeColor', 'none', ...
             'FaceColor', 'interp'); % Oberflaeche
view(-20, 60); % netter Blickwinkel
return

function C = thaller(Z, R)
[Nx, Ny] = size(Z);
% Farbtabelle gemaess Thaller
h = mod(-angle(Z) / (2*pi), 1);
t = pi - 2*atan2(abs(Z), R);
l = 1 - t / pi;

% HSV Tupel
C = ones(Nx, Ny, 3);
C(:,:,1) = h;
C(:,:,3) = l;
return