I've tried to solve the following PDE system with Neumann and Robin boundary conditions numerically for a long time, but it just doesn't work:
$\frac{\partial c_+}{\partial t} = -v \frac{\partial c_+}{\partial x} + D\frac{\partial^2 c_+}{\partial x^2}$
$\frac{\partial c_-}{\partial t} = v \frac{\partial c_-}{\partial x} + D\frac{\partial^2 c_-}{\partial x^2}$
The boundary conditions:
$D\frac{\partial c_+}{\partial x}(L,t) = 0$; $D\frac{\partial c_-}{\partial x}(0,t) = 0$
$vc_+(0,t) - D\frac{\partial c_+}{\partial x}(0,t) = v(I_E-W_If(u_I(t)))$
$-vc_-(L,t) - D\frac{\partial c_-}{\partial x}(L,t) = -vW_Ef(u_E(t))$
with $f(u) = \frac{u^4}{2^4+u^4}$.
$u_I(t)$ and $u_E(t)$ are given by:
$\frac{du_E}{dt} = -\gamma u_E + \kappa_E c_+(L,t)$
$\frac{du_I}{dt} = -\gamma u_E + \kappa_I c_+(0,t)$
With given $u_I(t)$ and $u_E(t)$, I find a numerical solution of $c_+(x,t)$ and $c_-(x,t)$ with an implicit scheme. That works well.
The problem is the numerical solution of $u_I(t)$ and $u_E(t)$. Does anyone have an idea how to find the solution of $u_I(t)$ and $u_E(t)$ numerically?
