(defun xfga-wavefunction (x0 p0 a-sub-n epsilon epsilon-derivative phi-sub-u)
  "Returns the XFGA wavefunction for the given parameters per Garraway 2000."

 #'(lambda (x)
      (let* ((epsilon-magnitude (abs epsilon))
             (sum-length (length a-sub-n))
             (prefactor (* (/ 1 (sqrt (* (sqrt pi) epsilon)))
                           (exp (* #C (0.0 1.0)
                                   (+ (/ (* 0.5 epsilon-derivative
                                            (expt (- x x0) 2)) epsilon)
                                      (* x p0)
                                      (* -0.5 x0 p0)
                                      (* -1 phi-sub-U)))))))
        (loop for n from 0 below sum-length
           collecting (* (expt (/ epsilon-magnitude (* (sqrt 2) epsilon)) n)
                         (hermite n (/ (- x x0) epsilon-magnitude))
                         (/ 1 (sqrt (factorial n)))) into coefficients
           finally (return (* prefactor
                              (reduce '+ (mapcar '* a-sub-n coefficients))))))))