function u = heat(k, t, x, init, bdry)
% solve the 1D heat equation on the rectangle described by
% vectors x and t with u(t(1), x) = init and Dirichlet boundary
% conditions u(t, x(1)) = bdry(1), u(t, x(end)) = bdry(2).

J = length(x);
N = length(t);
dx = mean(diff(x));
dt = mean(diff(t));
s = k*dt/dx^2;

u = zeros(N,J);
u(1,:) = init;

for n = 1:N-1
	u(n+1,2:J-1) = s*(u(n,3:J) + u(n,1:J-2)) + (1-2*s)*u(n,2:J-1);
	u(n+1,1) = bdry(1);
	u(n+1,J) = bdry(2);
end