# Stabilization of solution to one-dimensional system of PDE

I am trying to solve numerically next PDE system: $$\frac{\partial c}{\partial t}=\epsilon\frac{\partial}{\partial x}(\frac{\partial c}{\partial x}+\rho\frac{\partial \varphi}{\partial x}+\frac{vc}{1-vc}\frac{\partial c}{\partial x})$$ $$\frac{\partial \rho}{\partial t}=\epsilon\frac{\partial}{\partial x}(\frac{\partial \rho}{\partial x}+c\frac{\partial \varphi}{\partial x}+\frac{v\rho}{1-vc}\frac{\partial c}{\partial x})$$ $$-\epsilon^2\frac{\partial \varphi^2}{\partial x^2}=\rho$$ with completely blocking boundary conditions at x= ±1; $v$ and $\epsilon$ are constants. However, at some point of calculation the oscillations near the boundaries occur due to the too big gradients near the boundaries. I have found some information that there are several techniques for handling numerical instabilities. (Petrov-Galerkin (SUPG) or Galerkin least-squares (GLS) method). Unfortunately all of them usually described for the case of simple Convection-Diusion equation. Can anybody help me how to adjust some of those methods (or maybe another method) for my case?

-
You can also try to ask this at scicomp.stackexchange.com –  András Bátkai Jun 19 '13 at 7:34
Could it happen that $1-vc$ vanish at some point? –  Andrew Jun 19 '13 at 16:41
Definitely such situation can happen, however it is not a problem. Not going into technical details I can say, that this system was numerically solved for a range of values (which corresponds to a small concentration ​​and voltage values in corresponding physical model) despite the fact that $(1−vc)$ tended to zero. Problems occur with increasing of values ​​of certain parameters. Experimental results indicate that the problem is a loss of stability. –  user35096 Jun 20 '13 at 5:16