function [U,X,T]=heatImplEuler(T,f,v,g0,g1,M,MT)

% [U,X,T]=heatImplEuler(T,f,v,g0,g1,M,MT)
% 
% Implicit Euler method for
%     u_t - u_xx = f(x,t)     0<x<1 0<t<T
%     u(x,0)=v(x)
%     u(0,t)=g0(t)  u(1,t)=g1(t)
% 
% Input:
%    T     : final time
%    f     : source term
%    v     : initial condition
%    g0,g1 : BCs
%    M     : number of intervals in x
%    MT    : number of intervals in t
%  Output:
%    U     : (M+1) x (MT+1) matrix with the solution
%    X,T   : space -time grid (as matrices)  

k=T/MT;
h=1/M;
lambda=k/h^2;
x=h*(1:M-1)'; 

I=[1:M-1 1:M-2 2:M-1];
J=[1:M-1 2:M-1 1:M-2];
D=[2*ones(1,M-1) -ones(1,M-2) -ones(1,M-2)];
matrix=speye(M-1)+lambda*sparse(I,J,D);

U=zeros(M-1,MT+1);       % Making room for the solution
U(:,1)=v(x);             % I.C.

for n=1:MT
    t=n*k;
    rhs=U(:,n)+k*f(x,t); 
    rhs([1 M-1])=rhs([1 M-1])+lambda*[g0(t); g1(t)];
    U(:,n+1)=matrix\rhs;
end

t=k*(0:MT); 
x=h*(0:M);
U=[g0(t);...     % These lines add the BCs (when no comp., we take BC over IC) 
   U;...
   g1(t)];
[T,X]=meshgrid(t,x);

return