Debugging

Run your routines with the main program in Appendix B. This program works with hydrogen. It first makes phi and Vxc for a zero charge density (testing that you have handled Rho[]=0 properly) and creates -1/r+phi+Vxc, which should be just -1/r because Rho=0. It then solves for the 1s state. From the density of the 1s state, it then computes Vxc and Depsxc, and finally $\int V_{xc} n \, dV$ and $\Delta E_{xc} = \int \Delta\epsilon_{xc} n \, dV$. You should find that your results converge to -0.3316096642 and 0.0772489584, respectively, for these quantities.