function u = heat_imp(u0,D,dt,n)
% function u = heat_imp(u0,D,dt,n)

nx = length(u0);
x = linspace(-1,1,nx);
dx = max(diff(x));

a = dt/(4*dx^2);
A = zeros(nx,nx);
A = A + (1 + 2*a)*eye(nx);
for i = 2:nx-1
   A(i,i-1) = -a;
   A(i,i+1) = -a;
end
A(1,1) = 1;
A(nx,nx) = 1;

A1 = inv(A);
u(:,1) = u0';

for i = 1:n
   u(:,i+1) = A1*u(:,i);
end

return