# Discretizing a parabolic PDE with finite volume method

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?