• Equation 16 should not have subscript on T's
• Also missing b_i in before \partial_i B
• We also ignored spatial gradient in 'n' in both terms (and in derivative of \chi_\perp)

This would give us two additional terms in the expansion of the first term in 'q'.

\left ( \partial_i n \right )\left ( \chi_\parallel - \chi_\perp \right) b_i b_j \partial_j T

-n\left (\partial_n\chi_\perp \right) \left ( \partial_i n \right )b_ib_j \partial_j T

But since \chi_\perp is linear in n, these two terms combine to give

\left ( \partial_i n \right )\left ( \chi_\parallel - 2\chi_\perp \right) b_i b_j \partial_j T

As well as modifying the last term to

\nabla \cdot \left ( n \chi_\perp \nabla T \right )

which can be expanded as

n \chi_\perp \partial_i \partial_i T + n \left ( \partial_B \chi_\perp \right ) \left (\partial_i B \right ) \left ( \partial_i T \right ) + 2 \chi_\perp \left ( \partial_i n \right ) \left (\partial_i T \right )

Then we can write

\partial_i T = \frac{T_{\hat{i}} - T_{-\hat{i}}}{2 \Delta x}

and

\partial_i \partial_j T = \frac{T_{\hat{i}+\hat{j}} - T_{\hat{i}-\hat{j}} - T_{-\hat{i}+\hat{j}} + T_{-\hat{i}-\hat{j}}}{4 \Delta x^2} \left ( 1-\delta_{ij} \right ) + \frac{T_{\hat{i}} - 2 T + T_{-\hat{i}}}{\Delta x^2} \delta_{ij}

and

\partial_i \partial_i T = \frac{T_{\hat{i}} - 2 T + T_{-\hat{i}}}{\Delta x^2}

And we need to linearize the product of derivatives

\left ( \partial_i T \right ) \left ( \partial_j T \right )

First we subsitute T = T' + \delta T and expand to first order in \delta T

\left ( \partial_i T' + \partial_i \delta T \right ) \left ( \partial_j T' + \partial_j \delta T \right ) = \left ( \partial_i T' \right ) \left (\partial_j T \right ) + \left ( \partial_j T' \right ) \left ( \partial_i T \right ) - \left ( \partial_i T' \right ) \left ( \partial_j T' \right )

Now if we group terms in the expansion

A = B_i \partial_i T + C_{ij} \partial_i \partial_j T + D_{ij} \left ( \partial_i T \right ) \left ( \partial_j T \right ) + E \partial_i \partial_i T

we can work out coefficients for the various stencil terms (note D_{ij}+D_{ji} = 2D_{ij} since D_{ij} is symmetric)

\alpha_0 = -\frac{1}{\Delta t} -\frac{2C_{ii}}{\Delta x^2} - \frac{2\delta_{ii}E}{\Delta x^2}

\alpha_{\pm \hat{i}} = \pm\frac{B_i}{2\Delta x} + \frac{C_{ii}}{\Delta x^2} \pm \frac{2D_{ij} \left ( T'_{\hat{j}} - T'_{-\hat{j}} \right )}{4 \Delta x^2} + \frac{E}{\Delta x^2}

\alpha_{\pm \hat{i} \pm \hat{j}} = \pm \pm \frac{C_{ij}}{4\Delta x^2}\left ( 1 - \delta_{ij} \right )

and the right hand side

A = -\frac{T'}{\delta t} + \frac{D_{ij}}{4 \Delta x^2}\left ( \partial_i T' \right ) \left ( \partial_j T' \right )