% 2(i) Thin solid film growing with slope selection
% PDE: u_t = -e^2*(u*u_xxx)_x-(u^3*u_x)_x
%  IC: u(x,0) = exp(-10x^2)
%      x in [-pi, pi] and e = 0.2

% Grid, variable coefficient, and initial data:
  N = 128;                          % number of grid points
  h = 2*pi/N;                       % spatial step
  x = h*(-N/2+1:N/2)';                   % spatial grid
  t = 0;                            % initial time value
  dt = 1/N;                         % stable time step
  k = [0:N/2-1 0 -N/2+1:-1];
  e = 5;
  
  % initial condition: define u = u(x,0) and uold = u(x,-dt)
  u = exp(-x.^2)';
  uold = exp(-(x-dt).^2)';

  
% Time-stepping by Crank-Nicolson
  tmax = 5;                     % end time value
  tplot = .05;                  % time step for plots
  clf, drawnow;                 % clear figure
  set(gcf,'renderer','zbuffer');% set rendering for fast plotting
  H = gcf;                      % get handle to current figure
  plotgap = round(tplot/dt);    % time gap between plots
  dt = tplot/plotgap;           % adjust dt to line up with time plots
  nplots = round(tmax/tplot);   % number of time values to plot
  data = [u; zeros(nplots,N)];  % approximate solutions for plotting
  tdata = t;                    % keep track of time values used
  h = waitbar(0,'please wait...'); % create waitbar
  counter = 0;
  for i = 1:nplots              % loop through all time steps to be plotted
    waitbar(i/nplots)           % update wait bar
    
    for n = 1:plotgap           % loop through all time steps between plots
      t = t+dt;                 % increment time
      counter = counter + 1;                        
      u_k = fft(u);
      uold_k =fft(uold);
      u_x = real(ifft(i*k.*u_k));
      uold_x = real(ifft(i*k.*uold_k));
      f_k = i*k.*(fft(u_x.^3-u_x));
      fold_k = i*k.*(fft(uold_x.^3-uold_x));
      fhalf = (3*f_k - fold_k)/2;
      l1 = (2/dt^2+k.^2/dt-1/2);
      l2 = (1/dt^2+k.^2/dt+1/2);
      unew = (fhalf + l1.*u_k- uold_k/dt^2)./l2;
      unew = real(ifft(unew));
      uold = u;
      u = unew;
  end

    data(i+1,:) = u;            % include solution for plotting
    tdata = [tdata; t];         % time of current solution
    if max(abs(u)>100)          % check for instability
        data = data(1:i+1,:);   % stop when instability dominates too much
        break 
    end
  end
  
  % Trefethen's plotting code
  waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
  axis([-pi pi 0 tmax min(min(data)).5]), xlabel x, ylabel t, zlabel u;
  grid off, close(h),figure(H);