function u = heatvc(k, t, x, init, bdry)
% solve the 1D heat equation with variable coefficient k 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(2:J-1).*(u(n,3:J) + u(n,1:J-2)) + ...
		(1 - 2*s(2:J-1)).*u(n,2:J-1) + ...
		0.25*(s(3:J) - s(1:J-2)).*(u(n,3:J) - u(n,1:J-2));
	u(n+1,1) = bdry(1);
	u(n+1,J) = bdry(2);
end