function u = heatdiscrete(k, x, t, init, bdry)
% solve the 1D heat equation on the rectangle described by
% vectors x and t with u(x, t(1)) = init and Dirichlet boundary
% conditions u(x(1), t) = bdry(1), u(x(end), t) = 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

