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