Residual vector $\vec g_{it+1}$:

Once the trial Rho[it+1] is computed, you must also compute the corresponding g[it+1] with a call such as

    getg(g[it+1],Rho[it+1],Z,lmax,nmax,nmaxmax,F,r,dr,N);