function [x,t,U] = solveHeatImpl(a,f,T,u0,g0,g1,h,tau)
% function [x,t,U] = solveHeatImpl(a,f,T,u0,g0,g1,h,tau)
% Solving the heat equation on rectangle [0,1]x[0,T] using Finite Differences Method - Implicit scheme
% The equation:
%
%    du/dt = a d^2u/dx^2 + f(x,t)
%
% a     - constant (properties of the material)
% f    	- annonymous function of variable x and t
% T     - upper bound for time
% u0    - initial condition (annonymous function of variable x)
% g0,g1	- boundary conditons (annonymous functions of variable t)
% h     - step in x
% tau   - step in time
% x,t  - knots in both directions including boundaries
% U    - approximate solution in the knots - matrix
%

N=round(1/h);
h=1/N;              % recounting of h
x=linspace(0,1,N+1);
xx=x(2:end-1);		% inner points
t=0;
U=u0(x);            % initial condition
UU=U(2:end-1);      % initial condition in inner points
UU=UU';
r=a*tau/h^2;
tt=t;

while tt<T
    tt=tt+tau; t=[t,tt];
    dB=1+2*r;lB=-r;uB=-r;
    fk=tau*f(xx,tt); fk=fk(:);
    B0=g0(tt);
    B1=g1(tt);
    fk(1)=fk(1)+r*B0;
    fk(end)=fk(end)+r*B1;
    UU=solve3diag(dB,lB,uB,UU+fk);
    U1=[B0 UU' B1];
    U=[U;U1];
end

end     % That's all, Folks!

