% Usage: [u x t] = heatimp(f,T,n,nt) or u = heatimp(f,T,n,nt)
% Heat equation on [0,1] with an implicit method
%
% Input:
% f - Matlab inline function f(x) for the initial condition
% T - final time moment
% n - number of subintervals (panels)
% nt - number of time intervals
%
% Output:
% u - computed solution
% x - nodes of the spacial grid
% t - time steps
%
% Examples:
% [u x t] = heatimp(@bump,1,100,100);          here 'bump' is a user-defined function in M-file
% [u x t] = heatimp(@bump,1,100,80);
% [u x t] = heatimp(@bump,1,100,99);

function [u x t] = heatimp(f,T,n,nt);
x = (0:n)/n;
t = (0:nt)*T/nt;
u = zeros(n+1,nt+1);
h = 1/n;
tau = T/nt;
c = tau/h^2;
c1=1+2*c;
A=diag(c1*ones(n+1,1))-diag(c*ones(n,1),1)-diag(c*ones(n,1),-1);
A(1,1)=1;
A(1,2)=0;
A(n+1,n+1)=1;
A(n+1,n)=0;
for i=1:n+1
    u(i,1) = f((i-1)*h);
end;
for k=1:nt-1
u(:,k+1) = A\u(:,k);
end;