% Solves u_t + c(x) u_x = 0, in [0,2*pi] with periodic boundary conditions
%
  N = 256; h = 2*pi / N; x = h*(1:N);
  dt = h / 4;
  c = 0.2 + sin(x-1).^2;

% Initial condition

  v = exp(-100*(x-1).^2);
  vold = exp(-100*(x-0.2*dt-1).^2);

% Leapfrog scheme

  t = 0;  tmax = 8;  tplot = 0.15; clf, drawnow, set(gcf, 'renderer', 'zbuffer')
  plotgap = round(tplot/dt); dt = tplot / plotgap;
  nplots = round(tmax/tplot);
  data = [v; zeros(nplots,N)]; tdata = t;
%
  for j=1:nplots
     for n=1:plotgap
	t = t + dt;
	v_hat = fft(v);
	w_hat = 1i * [0:N/2-1 0 -N/2+1:-1] .* v_hat;
	w = real(ifft(w_hat));
	vnew = vold - 2*dt*c .* w; vold = v; v = vnew;
     end
     data(j+1,:) = v; tdata = [tdata; t];
  end
%
  waterfall(x, tdata, data), view(10,50), colormap(1e-6*[1 1 1]);
  axis([0 2*pi 0 tmax 0 4]), ylabel 't', zlabel 'u', grid off;
