% 1(h) Klein-Gordon
% PDE: phi_tt - phi_xx + V'(phi) , V(phi) = .5*phi^2+s*phi^4
%      phi(x,0) = sech((x-2)/sqrt(1-c^2)), phi_t(x,0) = 0
%      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 = 3.84/N;                     % stable time step
  dt = 3.85/N;                     % unstable time step
  s = 10;                          % sigma for V'(phi)
  c = .98;                         % parameter of initial condition
  v = sech((x-2)/sqrt(1-c^2));     % initial condition
  vold = v;                        % ~ u(x,-dt) since u_t(x,0)=0
  
% Time-stepping by Leap Frog formula:
  tmax = 20;                       % 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
          f = v + 4*s*v.^3;
          vnew = 2*v-vold+dt^2*(w-f); % 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(-20,25),colormap(1e-6*[1 1 1]);
  axis([-pi pi 0 tmax min(min(data)) 1]), pbaspect([1 1 .5]),
  xlabel x, ylabel t, zlabel u, grid off, close(h);figure(H);