Implicit method for two coupled PDEs

I have two equations (coupled), with the variables $T_1$ and $T_2$ and the constant $T_0$, which are (when written unitless, i.e. without prefactors): $$\partial_t T_1 = 1-T_1^3+T_1+\nabla\left(\frac{1}{T_2^{2.5}+T_2^{2.4}}\right)\nabla T_1+\frac{1}{T_2^{2.5}+T_2^{2.4}}\nabla^2T_1$$ and $$\partial_tT_2=\left(\frac{T_1}{1+T_1^2}\left(T_0-T_2\right)+\nabla\left(\frac{1}{1+T_2+T_2^2}\right)\nabla T_2+\frac{1}{1+T_2+T_2^2}\nabla^2T_2\right)$$ Now my usual approach is to divide such equations into a linear and a nonlinear term, as for the first equation: $$\partial_t T_1 = \underbrace{1-T_1^3+T_1}_{\text{Nonlinear}}+\underbrace{\nabla\left(\frac{1}{T_2^{2.5}+T_2^{2.4}}\right)\nabla T_1+\frac{1}{T_2^{2.5}+T_2^{2.4}}\nabla^2T_1}_{\text{Linear}}$$ and then splitting the linear term up into $T_{x,0}$ and $T_{x,1}$ for the current and next step, as for a reduced equation $$\begin{split} \partial_tT_1&=\nabla T_1\\ &\Rightarrow\\ \frac{T_{1,1}-T_{1,0}}{dt}&=\frac{\nabla T_{1,1}+\nabla T_{1,0}}{2}\\ \left(\frac{1}{dt}-\nabla\right)T_{1,1}&=\left(\frac{1}{dt}+\nabla\right)T_{1,0} \end{split}$$ which then can be solved using a matrix approach and -solver. In order to extend it to $T_2$, I can create a vector of both variables, such that $$\vec{A}\begin{pmatrix} T_{1,1}\\ T_{2,1} \end{pmatrix}=\vec{B}\begin{pmatrix} T_{1,0}\\ T_{2,0} \end{pmatrix}+\begin{pmatrix} \text{Nonlinear}_{T_1}\\ \text{Nonlinear}_{T_2} \end{pmatrix}$$ But now I am stuck how to do that for $$\nabla\left(\frac{1}{T_2^{2.5}+T_2^{2.4}}\right)\nabla T_1$$ Of course I can split it up to $$\begin{split} \nabla\left(\frac{1}{T_2^{2.5}+T_2^{2.4}}\right)&=\frac{\nabla\left(\frac{1}{T_{2,1}^{2.5}+T_{2,1}^{2.4}}\right)+\nabla\left(\frac{1}{T_{2,0}^{2.5}+T_{2,0}^{2.4}}\right)}{2} \end{split}$$ but then I do not know how to integrate it into the matrix equations, after I do not know how to create a matrix $\vec{C}$, such that $$\vec{C}\cdot T_2=\nabla\left(\frac{1}{T_2^{2.5}+T_2^{2.4}}\right)$$ Is there such a matrix? Or is my whole approach skewed, and should be rethought?

• Why do you call it a linear term if it depends nonlinearly on $T_2$? I think you can linearize this operator (and all the other nonlinear things) by taking there $T_2$ from the previous (already computed) time layer. Or you can linearize more accurately, introducing the Jacobian. Jun 24 '17 at 0:12
• Can you give an example for that? Jun 24 '17 at 5:43
• Would method of line be a possible option? I like it a lot because you only have to worry about the spacial discretisation and let the ODE solver take care of the time-domain. Aug 18 '17 at 17:43

Linearization of nonlinear PDEs is usually done by using the Frechet derivative of the system. First, let us define $T=\left(\begin{array}{c}T_1\\T_2\end{array}\right)$. Then, in your coupled system you can always replace $T_1=\left(\begin{array}{cc}1&0\end{array}\right)T$ and $T_2=\left(\begin{array}{cc}0&1\end{array}\right)T$ and write it as a coupled system:
$$\partial_tT= F(T)$$. Now, you are looking for the frechet derivative of $F$. $$F_T\Delta:=\lim_{\epsilon\rightarrow0}\frac{F(T+\epsilon\Delta)-F(T)}{\epsilon}$$.