I want to discretize the following parabolic PDE:
$$u_t = \nabla\cdot(\alpha(x)\nabla u)- \beta u\\ x\in\Omega \subset \mathbb{R}^2\\ \partial_n u = 0\\ u(t,0) = u_0(x)\ge 0, \alpha(x)>0$$
Given the Neumann boundary condition above, one can integrate this PDE:
\begin{alignat}{2} \frac{\mathrm d}{\mathrm dt}\int_\Omega \! u\,\mathrm dx&=\phantom{-}\int_\Omega \! \nabla\cdot (a(x)\nabla u)\,\mathrm dx \:\:&-\int_\Omega \! \beta(x) u \,\mathrm dx \\ &=\phantom{-}\int_{\partial\Omega}\! a(x)\nabla u\cdot n \,\mathrm ds&-\int_\Omega \! \beta(x) u \, \mathrm dx \\ &=\phantom{-}\int_{\partial\Omega} a(x)\partial_n u \,\mathrm ds&-\int_\Omega \! \beta(x) u \, \mathrm dx \\ &=-\int_\Omega \! \beta(x) u \, \mathrm dx \end{alignat}
The form of the discretized linear system I want to get is
$$\frac{d\vec{u}}{dt}=A\vec{u}-\vec{c}$$
Unfortunately, I don't think this is exactly what I'm obtaining, despite I do think I use a correct discretization technique. To this end, please consider my steps below.
I discretize the two integrals using the finite volume method:
$$\frac{d}{dt}\sum\limits_{j\in n_i} u_j l_{ij}=-\sum\limits_{j\in n_i} \beta_{ij} u_j l_{ij},$$
where $n_i$ is the set of neighboring cells to the cell $V_j\subset \Omega$ and where $l_{ij}$ is the length of the boundary of $V_j$ between $u_i$ and $u_j$.
Now one can discretize it in time to obtain:
$$\sum\limits_{j\in n_i} \frac{u_j^{n+1}-u_j^n}{\Delta t}l_{ij}=-\sum\limits_{j\in n_i} \beta_{ij} u_j^n l_{ij}$$
So we obtain the linear system $$L\frac{d\vec{u}}{dt} = -A\vec{u},$$
where $L$ is the matrix containing $l_{ij}$ and each entry of $A$ contains $-\beta_{ij}l_{ij}$.
Am I correct? I think I'm quite confused as I don't see where each term of the matrix should go, and where the vector $\vec{c}$ is supposed to be. Moreover, what does one do with the $l_{ij}$ on the LHS?