Use your elements of b[itp] from the lu solver to form the final value of Rho[it+1] according to the formula (1). (See the note in Section 3.3 for hints on how to do this.)
This completes the implementation of the modified Broyden algorithm! You now should be able to reach complete convergence within about ten Broyden steps!!!
Hint: For larger atoms, the modified Broyden approach sometimes induces NaN errors. This happens because some the coefficients b[itp] can become negative and the formation of Rho[it+1] sometimes results in very small, negative densities (on the order of -1e-12). To prevent these negative densities from causing NaNs in your exc() routines, it is best to simply take the absolute value after forming Rho[it+1] with a loop like
for (k=0; k<=N; k++)
Rho[it+1][k]=fabs(Rho[it+1][k]);