code:

u(x, t); 0 < x < 1
du / dt = d2u / dx2
u(0, t) = u(1, t) = 0;
u(x, 0) = x(1 - x);

code:

clear;

h = 0.01;
tau = 0.25 * h * h;
x0 = 0.0;
x1 = 1.0;
n = fix((x1 - x0) / h) + 1;

grid = [x0 : h : x1];
prevU = zeros(n, 1);
for i = 1 : n
	x = (i - 1) * h;
	prevU(i) = x * (1.0 - x);
end
u = zeros(n, 1);

figure;
hold on;
plot(grid, prevU);
iteration = 1;
while 1
	u(1) = 0;
	u(n) = 0;
	for i = 2 : n - 1
		laplace = (prevU(i - 1) - 2 * prevU(i) + prevU(i + 1)) / (h * h);
		u(i) = prevU(i) + tau * laplace;
	end
	
	error = 0;
	for i = 2 : n - 1
		dudt = (u(i) - prevU(i)) / tau;
		d2udx2 = (u(i - 1) - 2 * u(i) + u(i + 1)) / (h * h);
		error = max(error, abs(dudt - d2udx2));
	end

	if mod(iteration, 100) == 0
		plot(grid, u);
		pause(0.025);
	end
	
	prevU = u;
	iteration = iteration + 1;
end