According to the 12-step program, it’s time to solve a 1D diffusion equation, aka Heat equation of the form,

$\displaystyle \frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2},$

where ${\nu}$ is the viscosity. The process is very similar to solving a 1D convection equation, but this time we need to approximate a second derivative. An easy way to get an expression that approximates a second derivative is by using Taylor series:

$\displaystyle u(x+h,t)=u(x,t)+hu'(x,t)+\frac{h^2}{2!}u''(x,t)+\cdots,$

and

$\displaystyle u(x-h,t)=u(x,t)-hu'(x,t)+\frac{h^2}{2!}u''(x,t)+\cdots.$

Now add them up and truncate to get,

$\displaystyle u(x+h,t)+u(x-h,t)=2u(x,t)+h^2u''(x,t).$

Since we are intereseted in the second derivative, we rewrite the last expression as

$\displaystyle u''(x,t)=\frac{1}{h^2}\left[u(x-h,t)-2u(x,t)+u(x+h,t)\right].$

The only part of the code that needs to change is the finite-difference loop itself:

// Finite-difference loop:
for (int it=1; it<=nt-1; it++)
{
for (int k=0; k<=nx-1; k++)
{
un[k][it-1] = u[k][it-1];
}
// Keep values of u at the 'ends' constant:
u[0][it] = un[0][it-1];
u[nx-1][it] = un[nx-1][it-1];
for (int i=1; i<=nx-2; i++)
{
u[i][it] = un[i][it-1] + vis*dt/dx/dx*(un[i+1][it-1]-2*un[i][it-1]+un[i-1][it-1]);
}
}


Imagine that we have a rod of length 2m where the initial temperature distribution is the step-function as defined in the code. If no heat is added, we expect the rod to cool over time. The shown plot supports this.