Debugging

Check that your output states are orthonormal and diagonalize $\mu$ by executing

[Psi, epsilon]=getPsi(W); %# Run getPsi
Psi'*O(Psi) %# Should be the identity
Psi'*H(Psi) %# Should be diagonal
epsilon %# Should match the diagonal elements of previous matrix