function A=heat(n,k)
%
%
%
T=0.02;
Nsteps=round(T/k);
U0=inifun(n);
splt(n,U0);
Ui=U0(2:end-1,2:end-1);

h=1/n;
r=k*n*n;
e=ones(n-1,1);
B = spdiags([e -4*e e], -1:1, n-1, n-1);
C = spdiags([e 0*e e], -1:1, n-1, n-1);
I = speye(n-1);
A=kron(I,B)+kron(C,I);
M=speye((n-1)^2)-r/2*A;
N=speye((n-1)^2)+r/2*A;

perm = colamd(M); % column minimum degree reordering
[LL,UU,PP] = lu(M(:,perm));

u=reshape(Ui',[(n-1)^2,1]);
for j=1:Nsteps,
    hs=N*u;
    y=LL \ (PP*hs);
    u(perm,1) = UU \ y;
%    u=M\(N*u);
    Ui=reshape(u,[n-1,n-1])';
    U=zeros(n+1);
    U(2:end-1,2:end-1)=Ui;
    splt(n,U);
end



