// evolve_demo -- evolve a wavepacket in a potential with time, given
//                the energy eigenvectors and eigenvalues as provided
//                by eigs

// arguments -- a -- eigenvalues, eigenvectors, and x values as provided
//                   by ho_demo.r, anho_demo.r, etc.
//              x0 -- initial wavefunction position
//              delx -- initial wavefunction width
//              nt -- number of time steps
//              dt -- time step interval

// output -- a matrix in which the columns represent the absolute square
//           of the wavefunction at successive times

evolve_demo = function(a,x0,delx,nt,dt)
{
  if (nargs != 5) {
    return "Usage: evolve_demo(a,x0,delx,nt,dt)";
  }
  energies = a.eigvals'; // energies are just eigenvalues of the hamiltonian
  q = a.eigvecs';        // eigenvectors transform from x to energy space
  x = a.x;               // grid cell positions
  n = length(x);         // get size of domain
  psi0 = exp( -(x - x0).*(x - x0)/delx^2);  // calc initial wavefunction
  psisq = real(psi0.*conj(psi0));   // get probability distribution
  phi = q*psi0;          // transform to energy space

  xpgopen();
  xpgclear();
  for (it in 1:nt) {     // loop on time
    time = it*dt;
    phi = phi.*exp(-1i*energies*dt);   // evolve wave packet in energy space
    psi1 = q'*phi;                     // back to x space
    prob = real(psi1.*conj(psi1));
    psisq = [psisq, prob]; // matrix of probs at successive times
    xpgwindow(1,0,0,15,10);
    xpgaxes(1,a.x[1],-a.x[1],6,"x",0,1,6,"probability");
    xpgline(1,5,a.x,prob);
    sprintf(title,"Time = %f",time);
    xpglabel(1,.8*a.x[1],.7,title);
    xpgpause();
    xpgclear();
  }
  xpgclose();

  return psisq;
}