1D heat equation

https://en.wikipedia.org/wiki/Diffusion_equation

https://en.wikipedia.org/wiki/Heat_equation

INTRODUCTION

The diffusion equation is a parabolic partial differential equation. In physics, it describes the macroscopic behavior of many micro-particles in Brownian motion, resulting from the random movements and collisions of the particles. In mathematics, it is related to Markov processes, such as random walks, and applied in many other fields.

If the diffusion coefficient is constant then the diffusion equation is identical to the heat equation.

Unknown function u(x,t), x in [0,L], t in [0,T],
u(x,t) is the density of the diffusing material (isotropic diffusion)
or the temperature at the point x and time t,

heat equation u_t = D * u_xx, 

D is a positive coefficient called the diffusion coefficient 
or the thermal diffusivity of the medium.
u_t is time derivative, 
u_xx is double spatial derivative.
Initial condition u(x,0) = f(x).
Boundary condition u(0,t) = g1(t) and u(L,t) = g2(t).

#          j  u(0,t)=g1(t),   t increasing -->
#  x=0  +--+--+--+--
#       |  |  |  |
# u(x,0)+--+--o--+--
# =f(x) |  |  |  |
#       +--+--o--o--  u(i,j)
#       |  |  |  |
#     i +--+--o--+--
#       |  |  |  |
#  x=L  +--+--+--+--
#          j  u(L,t)=g2(t)

x = i*dx, i = 0..Nx (Nx+1 points), Nx*dx = L,
t = j*dt, j = 0..Nt (Nt+1 points), Nt*dt = T,

u_t = [ u(i,j+1) - u(i,j) ] / dt   # forward difference O(dt)
u_xx = [ u(i-1,j) -2*u(i,j) + u(i+1,j) ] / (dx*dx)   # central difference O(dx^2)

r = D*dt / (dx*dx), formula is stable if r <= 0.5

u(i,j+1) - u(i,j) = r * [ u(i-1,j) -2*u(i,j) + u(i+1,j) ]

u(i,j+1) = r*u(i-1,j) + (1-2*r)*u(i,j) + r*u(i+1,j)   # j+1 point in time

Initial condition, t=0 (j=0)
for i in range(Nx+1):
    u(i,0) = f(i*dx)

Boundary condition, x=0 (i=0) and x=L (i=Nx)
for j in range(1,Nt+1):
    u(0,j) = g1(j*dt)
    u(Nx,j) = g2(j*dt)

Final formula (inside mesh)
for j in range(Nt):
    for i in range(1,Nx):
        u(i,j+1) = r*u(i-1,j) + (1-2*r)*u(i,j) + r*u(i+1,j)

# Faster final formula using vectorization
for j in range(Nt):
    u[1:-1,j+1] = r*u[:-2,j] + (1-2*r)*u[1:-1,j] + r*u[2:,j]