W = 0.5;
alpha = 0.0001;
n = 1:1024*5;
L = 64*5;
N = length(n)

c0 = 30
xstep = L / N;
x = (n - 1)*L/N;
k = 2 * pi * (n - 1)/(L);

omega = c0 * k .* (1 - alpha * k .* k);
offset = 10
dh0 = exp(-(x-10).*(x-10)/(2*W^2));

%norming
area = 0;
for ii = 2 : length(x)
    area = area + (dh0(ii - 1) + dh0(ii)) / 2 * xstep;
end
C = 1 / area
dh0 = dh0 * C;

dhk = fft(dh0)

NT = 100;
tstep = 0.1;
for tt = 1 : NT
    time = (tt - 1) * tstep;
    for ii = 1 : 1 + N/2
        dhkt(ii) = dhk(ii) * exp(-j * omega(ii) * time);
    end
    for ii = N / 2 + 2 : N
        dhkt(ii) = dhk(ii) * exp(j * omega(N - ii + 2) * time);
    end
    dhx = ifft(dhkt);
    dhxarray(tt, :) = dhx;
end
[XX, YY] = meshgrid(x, 0:tstep:tstep*(NT-1));
size(XX)
figure(1)
surf(XX, YY, dhxarray);
ylabel('Time in s')
xlabel('distance from 0')
zlabel('the wave')
figure(2)
plot(x, dhxarray(1,:), x, dhxarray(10,:), x, dhxarray(40, :), x, dhxarray(100,:))
title('Wave contours at timesteps 0s, 10s, 40s, 100s')
xlabel('distance from 0')
ylabel('wave')
