First of all got to clarify that I have a biology background, so some things might seem obvious to the average math person. I have searched a lot about a solution to this and I only come across simple, constant boundary conditions which is not the case in my example. Also, there is no analytical solution possibility. I am thoroughly explaining the problem in purpose, from a biologist's point of view.
Background: I am trying to simulate a system of coupled ODEs and PDEs concerning the diffusion of two molecules between two different cells. For the sake of simplicity, let's just be concerned with the diffusion PDE of a single molecule from one cell to another.
Assume that the concentration of this molecule inside the cells and in the extracellular material is expressed as $u_a$, $u_b$ and $u_e$ respectively. Cells are assumed to be just points of production and they are separated by a distance $L$. Diffusion in this monodimensional distance is governed by a diffusion constant $\Theta$. There is also a rate of release/absorption from the inside of the cell to the outside, due to cell membrane separation, which is expressed by $\eta$.
The rate of change for quantities $u_a$ and $u_b$ is known, in terms of ordinary differential equations, and depends on the synthesis(if present), degradation, and absorption/release of this quantity to the extracellular material. Specifically:
$$\frac{du_a}{dt}(t) = K_Q s(t) + \eta \Big(u_e(t) - u_a(t) \Big) - \gamma_i u_a(t) \quad : \textrm{$a$ cell}$$ $$\frac{du_b}{dt}(t) = \eta \Big(u_e(t) - u_b(t) \Big) - \gamma_i u_b(t) \quad : \textrm{$b$ cell}$$
where $K_Q$ is a constant, $s(t)$ is just the concentration of another activating species (which itself is governed by another ODE), $\gamma_i$ the intracellular degradation constant and $\eta$ the cell wall diffusion constant ($\eta = 2\ min^{-1}$). Cell wall is assumed as point "boundary wall", that just allows molecules to pass with this $\eta$ coefficient, depending on their inside/outside differences. In cell $b$, there isn't any production of the molecule - just in/out flow and internal degradation.
Problem: I want to model the diffusion of this substance throughout the monodimensional domain (line) between the two cells and I have a problem applying the correct boundary conditions. I consider the boundary conditions as the points right outside the cells. Initial conditions are known (assumed to be zero everywhere: $u_a(0)=0$, $u_b(0)=0$ and $u_e(x_i,0) = 0$.).
As known, diffusion of a substance in one dimension is governed by the general equation:
$$\frac{\partial{u_e}}{\partial{t}} = \Theta \frac{\partial^2{u_e}}{\partial{x^2}} + g(x,t)$$
where $\Theta$ is the constant diffusion coefficient and $g(x,t)$ a term of addition or removal of concentration in each point. It differs for the boundary points.
If we discretize the line of length $L$ between the cells in $N$ points, then we get three cases: a set of $N-2$ ODEs for the intermediate points $\frac{du_e^i}{dt}$, $i \in \{1,...,N-1\}$, and two boundary conditions for points $i=0$ and $i=N$ (right outside the cells a and b. My problem is that I don't understand what form should these boundaries take. They should take into account the intracellular and extracellular concentration of the molecule ($u_a$, $u_b$) in order to calculate how much substance is diffused outwards or inwards the cells. I understand that I cannot use Dirichlet boundaries since this concentration is varying through time. I tried using Neumann boundaries (i.e. setting $\frac{du_e^0}{dt}(x_0,t)$ and $\frac{du_e^N}{dt}(x_N,t)$ instead of $u(x_0,t)$ and $u(x_N,t)$). Should I include diffusion to the extracellular material and degradation of the substance in these boundaries? If I do this, I obtain the following system:
$$\begin{cases} \begin{gather*} \frac{du_e^i}{dt} = \frac{\Theta}{\Delta x^2} \Big(u_e^{i+1}(t) - 2u_e^i(t) + u_e^{i-1}(t) \Big) - \gamma_e u_e^{i}(t) \textrm{ , for } i \in \{1,...,N-1\} \\[1.5em] \frac{du_e^0}{dt} = \eta \Big( u_a (t) - u_e^0(t) \Big) + \frac{\Theta}{\Delta x^2} \Big(2u_e^1(t) - 2u_e^0(t) \Big) - \gamma_e u_e^0(t) \\[1.5em] \frac{du_e^N}{dt} = \eta \Big( u_b (t) - u_e^N (t) \Big) + \frac{\Theta}{\Delta x^2} \Big(2u_e^{N-1}(t) - 2u_e^N(t) \Big) - \gamma_e u_e^N(t) \end{gather*} \end{cases}$$
which doesn't converge - concentrations in the extracellular material tend to $\pm \infty$.
The alternative would be to ommit the extracellular diffusion and degradation terms for $i=0$ and $i=N$, but would this be physically accurate?
For solving this system, I use a low-level Forward-Euler solver I coded into python. I always follow the necessary $dt \leq \frac{\Delta x^2}{2 \Theta}$ for the timestep.
Goals: Would like to get some insight as to why my extracellular concentrations do not converge: is it a problem of the finite differences method, even when following the restriction, or of the boundary conditions I'm using? I would prefer to solve the system using the simple methods I have already coded, rather than trying more elaborate techniques (implicit solution or Crank-Nicolson scheme).
P.S. My discretization is based on this life-saving manual: http://hplgit.github.io/prog4comp/doc/pub/p4c-sphinx-Python/._pylight006.html
EDIT: for clarification, I have added the exact equations for cells $a$ and $b$ concerning the quantity $u$. The question then boils down to how to connect extracellular $u$ with the boundary point $u$? For example, at point $i=N$, the diffused $u_e$ all the way from cell $a$ needs to be passed on to $\frac{du_e^N}{dt}$ and this would involve the inclusion of a term considering $u_e^{N-1}$ - if not, the rate of change will always be zero and the amount constant. In what way to include this concentration/point is what I cannot find out.
EDIT2 Adding the full set of equations after request. Diffused quantities are $Q_j^1$ and $Q_j^2$ molecules (correspond to $u$), where index $j$ defines the locality of these quantities: $c$ for controller cell, $p$ for producer cell and $e$ for extracellular material. All values in $[x]$ refer to concentrations of species $x$, and correspond to functions of time: $[x](t)$. $[Ref](t)$ is assumed constant.
\begin{gather} \textrm{Cell C}=\begin{cases} \begin{split} \frac{d[A\!:\!Q^2]}{dt} & = \Big( \chi_{A:Q_2,r,0} + \chi_{A:Q_2,r} \frac{K_{r}^{n_r}}{K_{r}^{n_r} + [Ref]^{n_r}} \Big) \cdot \\ & \Big( \chi_{A:Q_2,a,0} + \chi_{A:Q_2,a} \frac{[Q_c^2]^{n_q}}{K_{q}^{n_q} + [Q_c^2]^{n_q}} \Big) - \gamma_{A:Q_2}[A\!:\!Q^2] \end{split} \\[3ex] \frac{d[B]}{dt} = \chi_{B_0} + \chi_B \frac{[A\!:\!Q^2]^{n_b}}{K_b^{n_b} + [A\!:\!Q^2]^{n_b}} - \gamma_B [B] \end{cases}\\ \textrm{Cell P}=\begin{cases} \frac{d[C]}{dt} = \chi_{C_0} + \chi_C \frac{[Q_p^1]^{n_c}}{K_c^{n_c} + [Q_p^1]^{n_c}} - \gamma_C [C] \\ \frac{d[D]}{dt} = \chi_{D_0} + \chi_D \frac{K_d^{n_d}}{K_d^{n_d}+[C]^{n_d}} - \gamma_D [D] \end{cases} \\ \frac{d[Q_c^1]}{dt} = K_{Q^1} [B] + \eta \Big( [Q_e^1] - [Q_c^1] \Big) - \gamma_i [Q_c^1]\\ \frac{d[Q_p^1]}{dt} = \eta \Big( [Q_e^1] - [Q_p^1] \Big) - \gamma_i [Q_p^1] \\ \frac{\partial [Q_e^1]}{\partial t} = \eta \Big( [Q_c^1] - [Q_e^1] \Big) + \eta \Big( [Q_p^1] - [Q_e^1] \Big) - \gamma_e [Q_e^1] + \Theta \nabla^2[Q_e^1] \\ \frac{d[Q_c^2]}{dt} = \eta \Big([Q_e^2] - [Q_c^2] \Big) - \gamma_i [Q_c^2] \\ \frac{d[Q_p^2]}{dt} = K_{Q^2} [D] + \eta \Big([Q_e^2] - [Q_p^2] \Big) - \gamma_i [Q_p^2] \\ \frac{\partial [Q_e^2]}{\partial t} = \eta \Big([Q_c^2] - [Q_e^2] \Big) + \eta \Big([Q_p^2] - [Q_e^2] \Big) - \gamma_e [Q_e^2] + \Theta \nabla^2[Q_e^2] \\ \end{gather}
Parameter values and more info concerning the model : https://pubs.acs.org/doi/suppl/10.1021/acssynbio.6b00220/suppl_file/sb6b00220_si_001.pdf
My previous answer was incorrect; here's a more careful derivation.
In the inter-cellular space we have as a starting point the law of diffusion $$J = -\Theta \frac{\partial u_e}{\partial x}$$ for flux $J$, from which it follows that $$\frac{\partial u_e}{\partial t} = \Theta \Delta u_e - \gamma_e u_e$$ on $0 < x < L$ and \begin{align*} - \Theta \frac{\partial u_e}{\partial x}(0) &= g(0)\\ \Theta \frac{\partial u_e}{\partial x}(L) &= g(L) \end{align*} at the boundary. Here $g(x)$ is the rate of concentration injection at the boundary, due to diffusion through the cell wall: $$g(0) = \eta \left[u_a - u_e(0)\right]$$ $$g(L) = \eta \left[u_b - u_e(L)\right]$$
The boundary condition on the left can then be discretized as $$\Theta \frac{u_e^1(t+\Delta t) - u_e^0(t + \Delta t)}{\Delta x} = - \eta\left[u_a(t) - u_e^0(t)\right]$$ and likewise for the right boundary.
Here's a C++ implementation, and an animation of what the concentration looks like over time (left-most and right-most points are $u_a$ and $u_b$), for some arbitrarily-chosen values of the physical constants. Note that you need a quite small time step for the simulation to be stable. Switching to an implicit discretization will dramatically improve stability.