% 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 = .74/N;       % time step

  c = .2 + sin(x+2).^2;             % c(x)
  v = exp(-100*(x+2).^2);           % ~ u(x,0)
  %vold = exp(-100*(x-.2*dt+2).^2);  % ~ u(x,-dt)
  % For 2nd-order Adams-Bashforth method
  v_hat = fft(v);             % differntiation via FFT
  w_hat = 1i*[0:N/2-1 0 -N/2+1:-1] .* vold_hat; % discrete diff. of vold
  w = real(ifft(w_hat)); 
  vold = v + dt*c.*wold;
  
% Time-stepping by 2nd-order Adams-Bashford formula:
  tmax = 8;                    % end time value
  tplot = .15;                 % time step for plots
  clf, drawnow, set(gcf,'renderer','zbuffer') % graphics rendering info
  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
      % 2nd-order Adams-Bashford approximation
      v = vold - dt/2*c.*(3*w - wold); 
      vold = v; 
      wold = w;
    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
  
  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);