\rho c_v \frac{\partial T}{\partial t} = \nabla \left ( \kappa T^{n} \nabla T \right )

\rho c_v \frac{\partial T}{\partial t} =\nabla^2 \left ( \frac{1}{n+1} \kappa T^{n+1} \right )

For backward euler we evaluate the RHS at time j+1 although we still have to linearize

T_i^{{j+1}^{n+1}} = T_i^{j^{n+1}} + \left (n + 1 \right ) T_i^{j^n} \left ( T_i^{j+1} - T_i^j \right )

or equivalently

T_i^{{j+1}^{n+1}} = T_i^{j^n} \left (T_i^{j+1} + n \left (T_i^{j+1} - T_i^{j} \right ) \right )

begin{eqnarray} \rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & = & \frac{\kappa \Delta t}{\Delta x^2}\left (T_{i+1}^{j^n} T_{i+1}^{j+1} - 2 T_{i}^{j^n} T_{i}^{j+1} + T_{i-1}^{j^n} T_{i-1}^{j+1} \right ) \\ & - & \frac{n \kappa \Delta t}{ \left ( n+ 1 \right ) \Delta x^2} \left (T_{i+1}^{j^{n+1}} - 2 T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right ) \end{eqnarray}

Now to make the scheme 2nd order in time we can set T_i^{j+1}=\frac{T_i^{j}+T_i^{j+1}}{2} and we get

begin{eqnarray} \rho_i c_v \left ( T_i^{j+1}-T_i^j \right ) & = & \frac{\kappa \Delta t}{2 \Delta x^2}\left (T_{i+1}^{j^n} T_{i+1}^{j+1} - 2 T_{i}^{j^n} T_{i}^{j+1} + T_{i-1}^{j^n} T_{i-1}^{j+1} \right ) \\ & - & \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} \left (T_{i+1}^{j^{n+1}} - 2 T_{i}^{j^{n+1}} + T_{i-1}^{j^{n+1}} \right ) \end{eqnarray}

Where we can identify the left and right coefficients as

A_\pm = -\frac{\kappa \Delta t}{2 \Delta x^2}T_{i\pm1}^{j^n}

and the center coefficient as

A_0 = 2 \frac{\kappa \Delta t}{ 2\Delta x^2}T_{i\pm1}^{j^n} + \rho c_v

and the RHS left and right contributions

B_\pm = -\frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}}

and the RHS source term contribution as

B_0 = 2 \frac{\left (n - 1 \right ) \kappa \Delta t}{ 2 \left ( n+ 1 \right ) \Delta x^2} T_{i\pm 1}^{j^{n+1}} + \rho c_v T_i^j

and finally in the code the values for A and B are negated:

    kx = (0.5*kappa1*dt_diff/(dx**2))
    CALL getTemp(Info,i,j,k,T)
    stencil_fixed(0) = -2.0*REAL(nDim,8)*kx*T(0)**ndiff
    DO l = 1, 2*nDim
       stencil_fixed(l) = kx*T(l)**ndiff
    END DO  
    source = ((ndiff-1.0)/(ndiff+1.0))*dot_product(T(0:2*nDim),stencil_fixed(0:2*nDim)) - T(0)*gamma7*Info%q(i,j,k,1)
    stencil_fixed(0) = stencil_fixed(0) - gamma7*Info%q(i,j,k,1)

also note that the equation is multipled by the mean molecular weight… which is why c_v \rightarrow c_v \times \chi = \gamma_7 and \kappa \rightarrow \kappa \times \chi = \kappa_1