% 1(a) variable coefficient wave equation
% PDE: u_t + c(x)*u_x = 0 
%      c(x) = .2+sin(x+2).^2, u(x,0) = exp(-100*(x+2).^2)
%      x in [-pi, pi]

% 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.82/N;                      % stable time step
  dt = 1.83/N;                      % unstable time step

  c = .2 + sin(x+2).^2;             % c(x)
  v = exp(-100*(x+2).^2);           % ~ u(x,0)
  vold = exp(-100*(x-0*dt+2).^2);   % ~ u(x,-dt)
  
% Time-stepping by leap frog formula:
  tmax = 8;                     % end time value
  tplot = .15;                  % 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 = [v; zeros(nplots,N)];  % approximate solutions for plotting
  tdata = t;                    % keep track of time values used
  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 0 -N/2+1:-1] .* v_hat; % discrete diff. of v
      w = real(ifft(w_hat));    % ~ u_x
      vnew = vold - 2*dt*c.*w;  % leap frog approximation
      vold = v;                 % update v_{n-1}
      v = vnew;                 % update v_n
  end
    data(i+1,:) = v;            % include solution for plotting
    tdata = [tdata; t];         % time of current solution
    if max(abs(v)>2.5)          % 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 -3 3]), xlabel x, ylabel t, zlabel u, grid off
  close(h),figure(H);