Discretising a PDE for MATLAB using upwind scheme 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$$
 A: 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.
