%1D epitaxial growth
%    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    %% Fourier spectral methods for solving film epitaxy equation
%    %% h_t = (h_x^3 - h_x - h_{xxx})_x
%    %% with 12-periodic boundary condition
%    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
N=input('enter the number of points N:');
Tmax=input('the final time:');
L=6;
h=2*pi/N;
x=h*(1:N)';
x=L*(x-pi)/pi;
xx=-L:0.01:L;
t=0;
dt=0.002;
u=0.1*(sin(pi*x/2)+sin(2*pi*x/3)+sin(pi*x));
column1=[0 0.5*(-1).^(1:N-1).*cot((1:N-1)*h/2)]';
D1=pi/L*toeplitz(column1,column1([1 N:-1:2]));
column2=[-1/6-N^2/12 -0.5*(-1).^(1:N-1)./sin((1:N-1)*h/2).^2]';
D2=(pi/L)^2*toeplitz(column2,column2);
column3=[0 -N^2/8.*(-1).^(1:N-1)./tan((1:N-1)*h/2)+0.75*(-1).^(1:N-1).*cos((1:N-1)*h/2)./sin((1:N-1)*h/2).^3]';
D3=(pi/L)^3*toeplitz(column3,column3([1 N:-1:2]));
column4=[1/80*N^4-1/30+1/12*N^2 (N^2/4-1/2)*(-1).^(1:N-1)./sin((1:N-1)*h/2).^2-1.5*(-1).^(1:N-1).*cos((1:N-1)*h/2).^2./sin((1:N-1)*h/2).^4]';
D4=(pi/L)^4*toeplitz(column4,column4);
M=moviein(round(Tmax/100/dt));
counter=0;
%computation using RK
while t<Tmax
    k1=dt*(3*(D1*u).^2.*(D2*u)-D2*u-D4*u);
    k2=dt*(3*(D1*(u+k1/2)).^2.*(D2*(u+k1/2))-D2*(u+k1/2)-D4*(u+k1/2));
    k3=dt*(3*(D1*(u+k2/2)).^2.*(D2*(u+k2/2))-D2*(u+k2/2)-D4*(u+k2/2));
    k4=dt*(3*(D1*(u+k3)).^2.*(D2*(u+k3))-D2*(u+k3)-D4*(u+k3));
    u=u+(k1+2*k2+2*k3+k4)/6;
    t=t+dt;
%end
uu=interp1([-L; x],[u(N); u],xx,'spline');
  freq=20
    if mod(counter,freq)==0
        plot(xx,uu)
        text(4,max(abs(u))/2,['t=' num2str(t)])
        grid on
        M=[M getframe];
    end
    counter=counter+1;
end
movie(M)