import numpy as np import matplotlib.pyplot as plt

Parameters

nx, ny = 100, 100 # grid points dx = dy = 1.0 # spatial resolution dt = 0.1 # time step D = 1.0 # diffusion coefficient steps = 200 # number of time steps

Initial condition: a gaussian bump in the center

x = np.arange(nx) y = np.arange(ny) X, Y = np.meshgrid(x, y, indexing='ij') u = np.exp(-0.01*((X - nx/2)**2 + (Y - ny/2)**2))

def laplacian_2d(f): return (np.roll(f, -1, axis=0) + np.roll(f, 1, axis=0) + np.roll(f, -1, axis=1) + np.roll(f, 1, axis=1) - 4*f) / dx**2

for _ in range(steps): lap = laplacian_2d(u) u += dt * D * lap

Plot final state

plt.imshow(u, origin='lower') plt.colorbar() plt.title("2D diffusion after {} steps".format(steps)) plt.show()