Why negative carrier densities

Hi,

I’m trying to solve de DD model using n, p and V as main variables and FE method to solve the equations. I’m using Gummel’s iterative method to simulate a 2D PN junction but I’m getting negative densities at equlibrium in the first iteration even when the initial guest is pretty good. I think this is a common problem but i don’t understand why this happen and i would like to know if there is a “easy” way to avoid it.