clear all
global kvec filt epsilon
N=256;
K=2*floor(N/3);
kvec=[0:K/2,zeros(1,N-K-1),-K/2:-1];
filt=[ones(1,K/2+1),zeros(1,N-K-1),ones(1,K/2)];
epsilon=0.01;
Tfin=1;

x=[0:2*pi/N:2*pi-2*pi/N];
load('initial.dat');
y=initial(:,1)';
a0=ifft(y).*filt;
yold=y;
steps=4096;

h=Tfin/steps;
at=a0;
n=0;
for j=1:steps
    k1=h*f(at);
    k2=h*f(at+k1);
    at=at+(k1+k2)/2;
    n=n+1;
end
yt=real(fft(at));
plot(x,yt);
