Debugging

As a quick test, the operator H() should be Hermitian. This means that $(\vec a^\dagger H \vec b)^* = \vec b^\dagger H \vec a$ for any vectors $\vec a$ and $\vec b$. You may check this with

a=randn(prod(S),1)+i*randn(prod(S),1); %# Two random vectors
b=randn(prod(S),1)+i*randn(prod(S),1); %# Two random vectors
conj(a'*H(b))
b'*H(a)

Hint: Be sure to have run ``setup.m'' and ``sch.m'' recently so that all needed variables are defined.