$$ \frac {\partial C}{\partial t} = D \frac {\partial^2 C}{\partial x^2} - w \frac {\partial C}{\partial x}$$

Discretise the equation

$$\frac {C^{t+1}_{x} - C^{t}_{x}}{\Delta t} = \frac {D}{2} \Bigg(\frac {C^{t+1}_{x+1} - 2C^{t+1}_{x} + C^{t+1}_{x-1}}{\Delta x^2} + \frac {C^{t}_{x+1} - 2C^{t}_{x} + C^{t}_{x-1}}{\Delta x^2} \Bigg) -\frac {w}{2} \Bigg( \frac {C^{t+1}_{x+1} - C^{t+1}_{x-1}}{2 \Delta x} + \frac {C^{t}_{x+1} - C^{t}_{x-1}}{2 \Delta x} \Bigg)$$$$C^{t+1}_{x} = C^{t}_{x} + \frac {D \Delta t}{2 \Delta x^2} \Bigg( C^{t+1}_{x+1} - 2C^{t+1}_{x} + C^{t+1}_{x-1} + C^{t}_{x+1} - 2C^{t}_{x} + C^{t}_{x-1} \Bigg) - \frac {w \Delta t}{4\Delta x} \Bigg(C^{t+1}_{x+1} - C^{t+1}_{x-1} + C^{t}_{x+1} - C^{t}_{x-1} \Bigg)$$

if $\frac {w \Delta t}{4\Delta x} = \sigma$ and $\frac {D \Delta t}{2 \Delta x^2} = \lambda$

$$C^{t+1}_{x} = C^{t}_{x} + \lambda C^{t+1}_{x+1} - 2\lambda C^{t+1}_{x} + \lambda C^{t+1}_{x-1} + \lambda C^{t}_{x+1} - 2\lambda C^{t}_{x} + \lambda C^{t}_{x-1} - \sigma C^{t+1}_{x+1} + \sigma C^{t+1}_{x-1} - \sigma C^{t}_{x+1} + \sigma C^{t}_{x-1}$$

and rearranging,

$$-(\sigma + \lambda)C^{t+1}_{x-1} + (1 + 2\lambda)C^{t+1}_{x} + (\sigma - \lambda)C^{t+1}_{x+1} = (\sigma + \lambda) C^{t}_{x-1} + (1 - 2\lambda)C^{t}_{x} + (\lambda - \sigma) C^{t}_{x+1}$$

gives the discretised equation for $x > 1, x < L$.

$$-(\sigma + \lambda)C^{t+1}_{-1} + (1 + 2\lambda)C^{t+1}_{0} + (\sigma - \lambda)C^{t+1}_{1} = (\sigma + \lambda) C^{t}_{-1} + (1 - 2\lambda)C^{t}_{0} + (\lambda - \sigma) C^{t}_{1}$$$$-(\sigma + \lambda)C^{t+1}_{0} + (1 + 2\lambda)C^{t+1}_{0} + (\sigma - \lambda)C^{t+1}_{1} = (\sigma + \lambda)C^{t}_{0} + (1 - 2\lambda)C^{t}_{0} + (\lambda - \sigma) C^{t}_{1}$$$$(1 - \sigma + \lambda)C^{t+1}_{0} + (\sigma - \lambda)C^{t+1}_{1} = (1 + \sigma - \lambda)C^{t}_{0} + (\lambda - \sigma) C^{t}_{1}$$$$-(\sigma + \lambda)C^{t+1}_{L-2} + (1 + 2\lambda)C^{t+1}_{L-1} + (\sigma - \lambda)C^{t+1}_{L} = (\sigma + \lambda) C^{t}_{L-2} + (1 - 2\lambda)C^{t}_{L-1} + (\lambda - \sigma) C^{t}_{L}$$$$-(\sigma + \lambda)C^{t+1}_{L-2} + (1 + 2\lambda)C^{t+1}_{L-1} + (\sigma - \lambda)C^{t+1}_{L-1} = (\sigma + \lambda)C^{t}_{L-2} + (1 - 2\lambda)C^{t}_{L-1} + (\lambda - \sigma) C^{t}_{L-1}$$$$-(\sigma + \lambda)C^{t+1}_{L-2} + (1 + \sigma + \lambda)C^{t+1}_{L-1} = (\sigma + \lambda)C^{t}_{L-2} + (1 - \sigma - \lambda)C^{t}_{L-1}$$
$$ \begin{bmatrix} 1-\sigma + \lambda & \sigma - \lambda & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\\\ -(\sigma + \lambda) & 1 + 2 \lambda & \sigma - \lambda & 0 & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & -(\sigma + \lambda) & 1 + 2 \lambda & \sigma - \lambda & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & 0 & \ddots & \ddots & \ddots & \ddots & \ddots & 0 & 0 \\\\ 0 & 0 & 0 & 0 & \cdots & 0 & -(\sigma + \lambda) & 1 + 2 \lambda & \sigma - \lambda \\\\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & -(\sigma + \lambda) & 1+\sigma+\lambda \end{bmatrix} \begin{bmatrix} C_0^{t+1} \\\\ C_1^{t+1} \\\\ C_2^{t+1} \\\\ \vdots \\\\ C_{L-2}^{t+1} \\\\ C_{L-1}^{t+1} \end{bmatrix} = \begin{bmatrix} 1+\sigma-\lambda & \lambda-\sigma & 0 & 0 & \cdots & 0 & 0 & 0 & 0\\\\ \sigma + \lambda & 1-2\lambda & \lambda-\sigma & 0 & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & \sigma + \lambda & 1-2\lambda & \lambda-\sigma & \cdots & 0 & 0 & 0 & 0 \\\\ 0 & 0 & \ddots & \ddots & \ddots & \ddots & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 & 0 & \sigma + \lambda & 1-2\lambda & \lambda-\sigma & 0 \\\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma+\lambda & 1-\sigma-\lambda \end{bmatrix} \begin{bmatrix} C_0^{t} \\\\ C_1^{t} \\\\ C_2^{t} \\\\ \vdots \\\\ C_{L-2}^{t} \\\\ C_{L-1}^{t} \end{bmatrix} $$

In [ ]: