Hello,

First just a bit of background on my application. I am attempting to implement the Chebychev time propagator method using the FFT technique to evaluate the Hamiltonian (the d^2/dx^2 operator to be specific).

I decided to start simply with the ground stationary state of the SHO (ie a real gaussian in dimensionless units) and repeatedly operate by H, which should return the same state scaled by the eigenvalue (0.5 in this case), ie

do i=1,large
psi = H psi
psi=psi/.5
enddo

I found that my wavefunction quickly "blew up." Upon closer
inspection I saw that even on the first operation that although my eigenvalue was very good, the wavefunction itself at the 0 and N-1 grid points was the incorrect sign (although nearly zero). After the next iteration the 0,2 and N-1, N-3 points had the wrong sign and the magnitude of these points was growing quickly. Soon the entire wavefunction suffered from this problem.

I believe that the problem must be related to multiplying by (-omega^2) in the "frequency" space, but I am not certain what exactly I am doing wrong.

I am making use of the dfour1 routine. Another question I have is that my H(f) yields both real and imaginary portions, if anyone could tell me how fix this issue I would be quite happy!

Any comments/suggestions would be greatly appreciated! Thanks