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 
and 
. You should find that your results converge to
-0.3316096642 and 0.0772489584, respectively, for these quantities.