sigma = 1.0;
xmin = -10.0;
xmax = 10.0;
npts = 512;
nstates = 14;

dx = (xmax-xmin)/npts;
x = xmin + dx*(0:npts-1);

% -- initial state/wavefunction
psi_init = exp(-0.5*x.^2/sigma^2)*(pi*sigma^2)^(-0.25);

psi = zeros(nstates,npts);             % -- list to store oscillator states
psi(1,:) = psi_init;
for nn=2:nstates
     psi(nn,:) = raising_psi(psi(nn-1,:),xmin,xmax,npts,sigma);
end

function adag_fn = raising_psi(previous_fn,a,b,n,sigma)     % -- raising_operator
 dx = (b-a)/n;
 x = a + dx*(0:n-1);

 % -- going into fourier space using 'fft' library
 fwd_fft = fft(previous_fn);
 k = (2*pi/(b-a))*[0:n/2-1,0,-n/2+1:-1];
 dfk_dx = 1i*k.*fwd_fft;                     % -- 1st derivative in fourier space
 df_dx = ifft(dfk_dx);                       % -- back into position space

 % -- a^dagger acting on the previous state
 adag_fn = (x.*previous_fn/sigma - sigma*real(df_dx))/sqrt(2);
 norm_fn = adag_fn*transpose(adag_fn);
 adag_fn = adag_fn/sqrt(norm_fn);            % -- normalization
end