% p1f.m - Nonlinear optics
%         ieu_t + e^2u_xx + (1-|u|^2)|u|^2u = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(iek^2t)*u_hat.
%         u(x,0) = exp(-2x^2), e = 0.3

% Set up grid and two-soliton initial data:
  N = 128;                           % number of spatial grid points
  h = 2*pi/N;                        % spatial grid size
  x = (2*pi/N)*(-N/2+1:N/2)';        % spatial grid
  dt = 2/N;                          % stable time step
  dt = 3/N;                          % unstable time step
  u = exp(-2*(x.^2));                % initial condition
  e = 0.3;                           % parameter epsilon
  t = 0;                             % initial time
  v = fft(u);                        % Fourier transform of u (u_hat)
  k = [0:N/2-1 0 -N/2+1:-1]';        % permuted k to match freq. of fft
  ikv = e*1i*k.^2;                   % scaling factor for 2nd partial deriv.

% Solve PDE and plot results:
  tmax = 8;                          % end time
  tplot = .15;                        % time step for plots
  clf, drawnow;                      % clear figure 
  set(gcf,'renderer','zbuffer')      % for fast plotting
  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
  udata = [u zeros(N,nplots)];       % approximate solutions for plotting
  tdata = t;                         % keep track of time values used
  H = gcf;                           % get handle to current figure
  clf, drawnow;                      % clear figure
  h = waitbar(0,'please wait...');   % create waitbar
  E = exp(t*ikv);                    % integrating factor at t
  Edt2 = exp(dt*ikv/2);              % integrating factor at dt/2 
  Edt = exp(dt*ikv);                 % integrating factor at dt
  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;                      % set time for solution 
      g = 1i/e;                     % scaling factor of RHS 

      % 4th-order Runge-Kutta method
      F1 = dt*g.*fft(f1f(v,k));
      F2 = dt*g.*Edt2.*fft(f1f(v+F1/2,k));
      F3 = dt*g.*Edt2.*fft(f1f(v+F2/2,k));
      F4 = dt*g.*Edt.*fft(f1f(v+F3,k));
      v = Edt.*v + (Edt.*F1 + 2*Edt2.*(F2+F3) + F4)/6;
    end
     
    u = ifft(v);                    % compute solution for plotting
    udata(:,i+1) = u;               % keep track of plotted data
    tdata = [tdata t];              % keep track of plotted time
    if max(abs(u)>2.5)              % check for instability
       udata = udata(:,1:i+1);      % save data before instability
       break 
    end
  end
  
  % Trefethen's code for plotting
  waterfall(x,tdata,abs(udata'),angle(udata'));  %for complex values
  view(-20,25), xlabel x, ylabel t, zlabel u, axis([-pi pi 0 tmax 0 1]) 
  grid off, set(gca,'ztick',[-1 1]), close(h), pbaspect([1 1 .5]),figure(H);