Debugging

First, for basic debugging, delete the statement ``gbl_f=2;'' from ``dft.m'' and add the statement

gbl_f=2*ones(Ns,1); %# Constant fillings as a vector

just after the definition of Ns in ``dft.m''. Rerun ``dft.m'' and verify that your output (with perhaps changes in the fourteenth or fifteenth digits) is identical to before, at least through the finite-difference test.

Finally, for a non-trivial test, replace the above specification for gbl_f with

gbl_f=[2;2/3;2/3;2/3]; %# Variable filling test

Your code will now produce different output, but should pass the finite-difference test. Running to convergence, you should find a lower total energy (There are fewer electrons and less mutual repulsion.), but still a degeneracy pattern consistent with one s state and three p states. I find E=18.380, epsilon=[4.683; 6.339; 6.339; 6.339].