Need Help Discretizing this Heat Equation This might be a fairly easy question but since I haven't done much numerical PDE's before I don't really know what to do. I know that for a normal heat equation i.e.
$\large \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}$ (for $D$ arbitrary) we are able to discretize this equation to obtain an implicit numerical form:
$\Large \frac{u^{k+1}_{j} - u^{k}_{j}}{\Delta t} = D \frac{u^{k+1}_{j+1} - 2u^{k+1}_{j} + u^{k+1}_{j-1}}{\Delta x^2}$
which can then be rewritten as
$\large u^k_j = (1 + 2S)u^{k+1}_{j} - S (u^{k+1}_{j+1} + u^{k+1}_{j-1})$ where $S = \frac{D\Delta t}{\Delta x^{2}}$
and we can then use this to solve the numerical PDE through the use of a tridiagonal matrix and by solving a matrix system.
However, say for instance I now have a modified heat equation such as:
${\Large \frac{\partial u}{\partial t} =  D \frac{1}{ \frac{\partial^2 u}{\partial x^2}} - D}$
If I try and discretize this equation, I now get 
$\Large \frac{u^{k+1}_{j} - u^{k}_{j}}{\Delta t} =D \frac{1}{\frac{u^{k+1}_{j+1} - 2u^{k+1}_{j} + u^{k+1}_{j-1}}{\Delta x^2}}$
Following the same method as I wrote above, I now get:
$\Large u^k_j = u^{k+1}_j - D \frac{1}{\frac{u^{k+1}_{j+1} - 2u^{k+1}_{j} + u^{k+1}_{j-1}}{\Delta x^2}}$
The problem I'm having now however is that I'm not sure how to write this into a program since I'm not able to split up the $u$ terms in a neat way to allow the computer to work with the matrix system and there doesn't seem like there's a way to create a tridiagonal matrix. Did I write this the correct way?
Thanks in advance.
 A: That's quite a strange PDE. I'd rather not call it a heat equation. Do you really need an implicit discretization for it? Anyway, you've obtained a set of nonlinear algebraic equations for $u^{k+1}_j$ and may utilize some nonlinear solver, like Newton's method.
Note that the system of equations is still tridiagonal, that is
$$
\left\{
\begin{gathered}
&f_1(u_1, u_2) = u_1^{k}\\
&f_2(u_1, u_2, u_3) = u_2^{k}\\
&f_3(u_2, u_3, u_4) = u_3^{k}\\
&\vdots\\
&f_j(u_{j-1}, u_{j}, u_{j+1}) = u_j^k\\
&\vdots\\
&f_{n-1}(u_{n-2}, u_{n-1}, u_{n}) = u_{n-1}^{k}\\
&f_{n}(u_{n-1}, u_{n}) = u_{n}^{k}
\end{gathered}
\right.
$$
I've removed $k+1$ index for simplicity.
The Newton's method linearizes this system around some approximate solution $u_j^{(s)}$. To start with, one can take $u^{k}_j$ as starting $s = 0$ approximation $u_j^{(0)}$.
Let $u^{(s+1)}_j = u^{(s)}_j + \delta u_j$. Then linearized equation 
$$
f_j(u_{j-1}, u_j, u_{j+1}) = u^k_j
$$
becomes
$$
f_j(u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) 
+ \frac{\partial f_j}{\partial u_{j-1}} (u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) \delta u_{j-1}
+ \frac{\partial f_j}{\partial u_{j}} (u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) \delta u_{j}
+ \frac{\partial f_j}{\partial u_{j+1}} (u_{j-1}^{(s)}, u_j^{(s)}, u_{j+1}^{(s)}) \delta u_{j+1}
= u_j^k.
$$
Assembling that altogether we get a linear system for $\boldsymbol \delta \mathbf u$
$$
\mathbf f^{(s)} + J^{(s)} \boldsymbol \delta \mathbf u = \mathbf u^{k}
$$
with $J$ being a tridiagonal Jacobi matrix
$$
J = \frac{\partial \mathbf f}{\partial \mathbf u}.
$$
So the time stepping from $k$ to $k+1$ algorithm becomes:


*

*Use $\mathbf u^k$ as initial $s = 0$ approximation for $\mathbf u^{k+1}$

*Compute $\mathbf f^{(s)}$ Jacobi matrix $J^{(s)}$ using last know approximation

*Solve linear system for $\boldsymbol \delta \mathbf u$

*Use $\mathbf u^{(s+1)} = \mathbf u^{(s)} + \delta \mathbf u$ as new approximation and increase $s$. If $\boldsymbol \delta \mathbf u$ is sufficiently small, use $\mathbf u^{(s+1)}$ as $\mathbf u^{k+1}$, else continue from step 2.

