function u = heat(k, t, x, init, bdry)
% Solve the 1D heat equation on the rectangle described by
% vectors x and t with initial condition 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