Numerical issues with logarithm

I need to solve a Poisson-equation that’s coupled with nonlinear equations containing logarithms.
That code is unfortunately really unstable because I keep getting negative values during iterations that make the logarithm diverge.

I got it working by cutting the logarithmic function off using a conditional statement and Taylor expansion at that point – however, this solution seems not to work when i extend this problem to pure Neumann conditions. Is there any way to restrict the value range of a function to positive values? Or are there other possibilities to avoid this problem?