I have to solve, using some numerical method, several equations similar to: $$\frac{\partial\,}{\partial t}\rho(\mathbf{x}, t) = \Delta\left(-\rho^2(\mathbf{x}, t) + \rho^5(\mathbf{x}, t)+\Delta(\rho(\mathbf{x}, t))\right)$$ where $\Delta$ is the 2-dimensional Laplacian $\Delta = \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$, and $\rho(\mathbf{x}, t)$ is a function with the constraint $\rho(\mathbf{x}, t) \geq 0$ for all $\mathbf{x}$ and all $t$.
I've tried mostly using finite-difference methods, but I cannot figure how to correctly implement the constraint $\rho(\mathbf{x}, t) \geq 0$.
Do you have any suggestion and/or reference?
Thank you for your attention