function [xh,uh]=heattheta(xspan,tspan,nstep,theta,mu,...
               u0,g,f,varargin)
%HEATTHETA risolve l'equazione del calore con il
%  theta-metodo.
%  [XH,UH]=HEATTHETA(XSPAN,TSPAN,NSTEP,MU,U0,G,F,THETA)
%  risolve l'equazione del calore
%  D U/DT - MU D^2U/DX^2 = F nel dominio
%  (XSPAN(1),XSPAN(2))x(TSPAN(1),TSPAN(2)) utilizzando
%  il theta-metodo con condizione iniziale U(X,0)=U0(X)
%  e condizioni al bordo di Dirichlet U(X,T)=G(X,T) per
%  X=XSPAN(1) e X=XSPAN(2). MU e' una costante positiva
%  F=F(X,T), G=G(X,T) e U0=U0(X) sono inline function o
%  anonymous function o function definite in M-files.
%  NSTEP(1) e' il n.ro di intervalli in spazio
%  NSTEP(2) e' il n.ro di intervalli in tempo
%  XH contiene i nodi della discretizzazione
%  UH contiene la soluzione numerica al tempo TSPAN(2).
%  [XH,UH]=HEATTHETA(XSPAN,TSPAN,NSTEP,MU,U0,G,F,...
%  THETA,P1,P2,...) passa i parametri opzionali
%  P1,P2,...to alle funzioni U0,G,F.

%  numero di intervalli lungo la variabile X,
%  NSTEP(2) e' il numero di intervalli temporali.
h  = (xspan(2)-xspan(1))/nstep(1);
dt = (tspan(2)-tspan(1))/nstep(2);
N = nstep(1)+1;
e = ones(N,1);
D = spdiags([-e 2*e -e],[-1,0,1],N,N);
I = speye(N);
A = I+mu*dt*theta*D/h^2;
An = I-mu*dt*(1-theta)*D/h^2;
A(1,:) = 0; A(1,1) = 1;
A(N,:) = 0; A(N,N) = 1;
xh = (linspace(xspan(1),xspan(2),N))';
fn = f(xh,tspan(1),varargin{:});
un = u0(xh,varargin{:});
[L,U]=lu(A);
for t = tspan(1)+dt:dt:tspan(2)
    fn1 = f(xh,t,varargin{:});
    rhs = An*un+dt*(theta*fn1+(1-theta)*fn);
    temp = g([xspan(1),xspan(2)],t,varargin{:});
    rhs([1,N]) = temp;
    uh = L\rhs; uh = U\uh; fn = fn1; un = uh;
end
return
