% 1(b) 2nd order wave
% PDE: u_tt = u_xx 
%      u(x,0) = exp(-10x^2), u_t(x,0) = 0
%      x in [-pi, pi]

% Grid, variable coefficient, and initial data:
  N = 80;                          % 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 = 3.9/N;                      % stable time step
  dt = 4/N;                        % unstable time step

  v = exp(-10*(x.^2));             % ~ u(x,0)
  vold = v;                        % ~ u(x,-dt) since u_t(x,0)=0
  
% Time-stepping by Leap Frog formula:
  tmax = 8;                        % end time value
  tplot = .075;                    % time step for plots
  clf, drawnow, set(gcf,'renderer','zbuffer') % graphics rendering info
  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 = [v; zeros(nplots,N)];     % 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
  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
          v_hat = fft(v);          % differntiation via FFT
          w_hat = (1i*[0:N/2-1 N/2 -N/2+1:-1]).^2 .* v_hat; % scale v_hat
          w = real(ifft(w_hat));   % ~ u_xx
          vnew = 2*v-vold+dt^2*w;  % leap frog approximation
          vold = v;                % update v_{n-1}
          v = vnew;                % update v_n
      end
      data(i+1,:) = v;             % spatial data
      tdata = [tdata; t];          % time data
      if max(abs(v)>2.5)           % check for instability
          data = data(1:i+1,:);    % save data before instability
          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 -2 2]), xlabel x, ylabel t, zlabel u, grid off
  close(h);figure(H);