I am trying to discretise the PDE: $$\phi \frac{\partial c}{\partial t}+\frac{\partial j}{\partial x}=0$$ where $c$ is a function of $x$ and $t$, and $j=qc-D(q)\frac{\partial c}{\partial x}$, $q$ is the Darcy flux (function of $x$ and $t$) and $D(q)$ is the coefficient of diffusion/dispersion, and $\phi$ is a constant.
over the interval $(x_{i-1/2}, x_{i+1/2})$
My Progress
$$\begin{align*} &\phi \frac{\partial c}{\partial t}+\frac{\partial j}{\partial x}=0 \\ &\Rightarrow \frac{\partial c}{\partial t}+\frac{\partial}{\partial x}\left(\frac{qc}{\phi}-D(q)\frac{\partial c}{\partial x}\right)=0 \\&\Rightarrow \frac{\partial}{\partial t}\int_{x_{i-1/2}}^{x_{i+1/2}}c \ dx \ +\left|\frac{qc}{\phi}-D(q)\frac{\partial c}{\partial x}\right|_{x_{i-1/2}}^{x_{i+1/2}}=0 \end{align*}$$
What can I do to simplify the expression in terms of the values of the $c_i, q_i$ terms? (I am trying to get matlab to run this iteratively).
Conventions
$x_i$ is the midpoint of the $i$th interval.
$\delta x$ is the length of the interval.
$$f_i=\frac1{\delta x}\int_{x_{i-1/2}}^{x_{i+1/2}}f \ dx$$
The PDE $$ \phi\partial_t c + \partial_x(qc) = \partial_x(D(q)\partial_x c) $$ is a linear advection-diffusion equation, with non-constant advection speed $q/\phi$ and non-constant diffusivity $D(q)/\phi$.
Now, let us integrate the last equation of the OP over $t\in [t_n , t_{n+1}]$ with $t_n= n\Delta t$. We have $$ c_{i}^{n+1} = c_i^n - \frac{\Delta t}{\Delta x} \left(F_{i+ 1/2}^n - F_{i- 1/2}^n\right) , \tag{*} $$ where $c_i^n = \frac{1}{\Delta x}\int_{x_{i-1/2}}^{x_{i+1/2}} c(x,t_n)\,\text{d}x$ denotes the cell averages at the time $t_n$, and $$ F_{i+ 1/2}^n = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \left.\left(\frac{qc - D(q) \partial_x c}{\phi}\right)\right|_{x = x_{i+ 1/2}}\, \text{d}t $$ denotes the flux at the cell's boundaries. The exact flux $F_{i+ 1/2}^n$ is then approximated by an appropriate numerical flux. Note that a $\phi$ in the denominator may be missing in the OP.
A straightforward explicit finite-difference approximation of the PDE based on central differences writes \begin{aligned} &\phi\frac{c_i^{n+1}-c_i^{n}}{\Delta t} + \frac{q_{i+1}^{n} c_{i+1}^{n}-q_{i-1}^{n} c_{i-1}^{n}}{2\,\Delta x} \\ &\quad = \frac{1}{\Delta x} \left( D(q_{i+1/2}^{n}) \frac{c_{i+1}^{n}-c_{i}^{n}}{\Delta x} - D(q_{i-1/2}^{n}) \frac{c_{i}^{n}-c_{i-1}^{n}}{\Delta x}\right) , \end{aligned} where $c_i^{n} \simeq c(i\Delta x,n\Delta t)$. This method is of the form $(*)$ with the numerical flux $$ F_{i+ 1/2}^n \simeq \frac{q_{i}^{n} c_{i}^{n}+q_{i+1}^{n} c_{i+1}^{n}}{2\phi} - D(q_{i+1/2}^{n}) \frac{c_{i+1}^{n}-c_{i}^{n}}{\phi\,\Delta x} \, . $$ A similar method based on upwind differences writes \begin{aligned} &\phi\frac{c_i^{n+1}-c_i^{n}}{\Delta t} + q_{i}^{n+}\frac{c_{i}^{n}-c_{i-1}^{n}}{\Delta x} + q_{i}^{n-}\frac{c_{i+1}^{n}-c_{i}^{n}}{\Delta x} \\ &\quad = \frac{1}{\Delta x} \left( D(q_{i+1/2}^{n}) \frac{c_{i+1}^{n}-c_{i}^{n}}{\Delta x} - D(q_{i-1/2}^{n}) \frac{c_{i}^{n}-c_{i-1}^{n}}{\Delta x}\right) , \end{aligned} where $q_{i}^{n+} = \max\lbrace q_{i}^{n},0\rbrace$ and $q_{i}^{n-} = \min\lbrace q_{i}^{n},0\rbrace$. Variants are possible, such as using $$ \frac{q_{i}^{n+}c_{i}^{n}-q_{i-1}^{n+}c_{i-1}^{n}}{\Delta x} + \frac{q_{i+1}^{n-}c_{i+1}^{n}-q_{i}^{n-}c_{i}^{n}}{\Delta x} $$ or even $$ q_{i-1/2}^{n+}\frac{c_{i}^{n}-c_{i-1}^{n}}{\Delta x} + q_{i+1/2}^{n-}\frac{c_{i+1}^{n}-c_{i}^{n}}{\Delta x} $$ for the upwind differencing term.