My colleagues and I are trying to numerically integrate a physics-based PDE (phase field model of diffusive phase separation) of the form $$ \frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x} \left [ a(\phi) \frac{\partial \phi}{\partial x} + b(\phi) \frac{\partial^{3} \phi}{\partial x^{3}} \right ] \tag{1} $$ where $$ a(\phi) = (1 - \phi)/N + \phi - 2 \chi \phi (1-\phi), $$ $$ b(\phi) = \kappa \phi (1-\phi), $$ and $\chi$, $\kappa$ and $N$ are real positive numbers.
There are several challenges to numerically integrating this system (e.g. stiffness), and we have typically used a pseudo-spectral method for the spatial derivatives and a semi-implicit, variable step-size Euler stepping scheme for the time integration (more details can be found here and here).
However, some of our most persistent challenges relate to the fact that there are physical bounds that constrain $\phi \in [0, 1]$. We find that if the equation is integrated with sufficient accuracy (i.e. if we use a sufficiently small time step), this is not a problem. However, at large values of $N$ and $\chi$ it becomes very difficult to reach the necessary accuracy.
Aside from simply increasing the accuracy of our time integration (e.g. smaller time step, higher order method, exponential time differencing) are there methods that we can use to stabilize our integration and keep $\phi$ between 0 and 1?