# Coupled linear PDE equations (2nd and 1st order) - Numerical Method

I am trying to solve a coupled reaction-diffusion equations, using Crank Nicolson (implicit Finite Differences Method). I know how to solve them separately, but not simultaneously.

\begin{align} \partial_t C_f (x, t) &= \partial_{xx} C_f(x, t) - (k_3 + k_5) C_f(x, t) + k_4 C_b(x, t) \\ \partial_t C_b(x, t) &= k_3 C_f(x, t) - k_4 C_b(x, t) \end{align}

I transform it to Crank Nicholson form and I obtain:

\begin{align} r &:= D \frac{\Delta t}{\Delta x^2}\\ -r C_{f, i-1}^{n+1} + 2 (1 + r) C_{f, i}^{n+1} - rC_{f, i+1}^{n+1} &= r C_{f, i-1}^n + 2 (1 -r)C_{f, i}^n -r C_{f, i+1}^n + \Delta t \big(k_4 C_{b, i}^n - k_3 C_{f, i}^n \Big) \\ C_{b, i}^{n+1} &= C_{b, i}^n(1 - \Delta t k_4) + \Delta t k_3 C_{f, i}^n \end{align}

Is the equation for $$C_b(x,t)$$ correct? How can I couple them? I know how to solve them separately applying boundary conditions (Dirichlet both sides for $$C_f$$), and Neumann both sides for $$C_b$$.

Thank you very much!

I am not sure what the capital "D" is in your definition of $$r$$. You have discretized the second equation with Forward Euler, not Crank-Nicholson which would read $$\frac{C_{b, i}^{n+1} - C_{b, i}^{n}}{\Delta t } = \frac{1}{2} \bigg[ k_3 \Big( C_{f, i}^{n +1} + C_{f, i}^{n} \Big) - k_4 \Big( C_{b, i}^{n +1} + C_{b, i}^{n} \Big) \bigg]$$ Your first equations also does not look correct to me (where is $$k_5$$?). It should be \begin{align}\frac{C_{f, i}^{n+1} - C_{f, i}^{n}}{\Delta t } = \frac{1}{2} \Bigg[& \frac{C_{f,i-1}^{n + 1} - 2 C^{n + 1}_{f,i} + C^{n + 1}_{f, i -1}}{\Delta x^2} + \frac{C_{f,i-1}^{n} - 2 C^{n}_{f,i} + C^{n}_{f, i -1}}{\Delta x^2} \\ & - (k_3 + k_5) \Big( C_{f, i}^{n+1} + C_{f, i}^n \Big) + k_4 \Big( C_{b, i}^{n+1} + C_{b, i}^n \Big) \Bigg] \end{align} instead. You can couple the system by grouping the $$C_{b, i}^n, C_{f, i}^n$$ into one vector which is then updated $$n \to n+1$$ by solving the linear system defined by the discretizaions. You might need to alter the system of equations to respect the boundary conditions manually.