function for mapping (x,y,z) indices to a linear index
Tsparse = sparse([],[],[],Nx*Ny*Nz, Nx*Ny*Nz, 7*(NxNz-2)); % kinetic part of Hamiltonian matrix: (d^2/dx^2 + d^2/dy^2 + d^2/dz^2); NOTE: we'll have 7*(Nx; % define an inline helper function for mapping (x,y