To give some background on the problem, I'm trying to solve the following PDE via finite differences:
$$\frac{\partial\sigma}{\partial t}=\frac{\partial}{\partial x}\left[ \kappa(x,t)\left(\frac{\partial\sigma}{\partial x}+G\right)\right]$$
After expanding the equation it takes the form of a convection diffusion equation:
$$\frac{\partial\sigma}{\partial t}=\kappa(x,t)\frac{\partial^2\sigma}{\partial x^2}+\frac{\partial\kappa(x,t)}{\partial x}\frac{\partial\sigma}{\partial x}+G\frac{\partial\kappa(x,t)}{\partial x}$$
It should also be noted that $k(x,t)\geq 0$, $G>0$, and $\frac{\partial\kappa(x,t)}{\partial x}\leq0$. All the coefficients vary in time and space, and in some regions the PDE is convection-dominated, which is known to make the problem more difficult. Right now, however, I've managed to get results which make physical sense using an upwind semi-discrete scheme; essentially, I discretize the equation in space using central differences on the diffusive (second order) term and forward differencing on the convective (first order) term. I then use a stiff ODE solver in MATLAB to solve the resulting system of ODEs.
The equation is also subject to Neumann conditions $\sigma_x(0,t)=\sigma_x(L,t)=-G$, with $G>0$, and this is where I'm running into some trouble. I'm using central differences to incorporate the Neumann conditions into the scheme; if $x_1,x_2,\dots,x_n$ are the mesh nodes, I use $\sigma_x(0,t)\approx (\sigma(x_2,t)-\sigma(x_0,t))/2\Delta x$ and plug the expression for $\sigma(x_0,t)$ into the scheme centered at $x_1$. A similar approach is taken for the other boundary condition.
The solution behavior seems to be somewhat accurate on the interior nodes, but at the $x=0$ boundary, the function $\sigma(0,t)$ either oscillates, or increases monotonically when it should decrease monotonically, which is a big problem. It doesn't grow without bound, but the oscillations/monotonic increases make the solution unphysical.
I have two questions: (1) Upwinding is one approach used to limit spurious oscillations at the interior nodes, but are there any techniques to limit spurious oscillations at boundaries with Neumann conditions? (2) Is there a way to impose some kind of monotonicity or non-positivity restriction on the evolution of the solution at $x=0$?