function [x, phi] = initialconditions(n)
  
  dx = .25;
  x = (0:n-1)*dx - 4;
  
  alpha = 1.;
  k0 = -1.;
  
  im = %i;
  phi = sqrt(2*alpha*alpha/%pi)*exp(im*k0*x).*exp(-alpha*alpha*x^2);
  
endfunction

function [v] = potentialwell(x,vpot)
  
  v = zeros(x);
  v(x>3 & x<5) = vpot;
  
endfunction

function [newphi] = evolve(x, phi, dt, vpot)
  
  im = %i;
  
  dx = x(2)-x(1);
  nx = max(size(x));
  
  v = potentialwell(x, vpot);
  myphi = phi;
  
  phifft = fft(myphi,-1);
  
  //calculate k
  
  dk = 1/(nx*dx);
  k = (0:nx-1)*dk;
  k(nx/2+1:nx) = k(nx/2+1:nx)-k(nx);
  
  propagator = exp(im*k^2*dt);
  newphifft = propagator.*phifft;
  newphi = fft(newphifft,+1);
  
endfunction

function [allt,x,allphi] = runspectral(n,ndt,dt,vpot)
  
  [x,phi] = initialconditions(n);
  nx = max(size(x));
  
  allt = (0:ndt)*dt;
  allphi = zeros(ndt+1,nx);
  
  allphi(1,:) = phi;
  
  for i=2:ndt+1 do
    allphi(i,:) = evolve(x,allphi(i-1,:), dt,vpot);
  end

endfunction