clear all
global kvec epsilon c
N=64;
x=[0:pi/N:2*pi-pi/N];
for i=[1:size(x,2)]
  y(i)=u0(x(i));
end
a0=ifft(y);
a0(N+1)=0;
kvec=[0:N-1,-N:-1];
c=.4; epsilon=0.01;
h=0.01;
at=a0;
for t=[0:h:100]
    k1=h*f(at);
    at=at+(k1+h*f(at+k1))/2;
    yt=fft(at);
    plot(x,yt);
    drawnow
end
