First this is a fantastic piece of coding. I learned a lot just looking at the code an figuring out how it was running. Second, anyone downloading this code should immediately download Ref. 1 by J.S. Blakemore, and if you're really serious you should by the book "Semiconductor Statistics" by Blakemore as well.
In the charge neutraility equation I believe that the degeneracy factor should be 1/gD rather than just gD, agree?
Also, the ND should be ND[1-F(E)] rather than just ND*F(E), agree?
Finally, thank you very much Ernesto. Nice Work!!