function [x,u]=heattheta(xspan,tspan,nstep,theta,mu,...
               u0,g,f,varargin)
%HEATTHETA risolve l'equazione del calore con il
%  theta-metodo.
%  [X,U]=HEATTHETA(XSPAN,TSPAN,NSTEP,THETA,MU,U0,G,F)
%  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. NSTEP(1) e' il
%  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;
x = linspace(xspan(1),xspan(2),N);
x = x';
fn = feval(f,x,tspan(1),varargin{:});
un = feval(u0,x,varargin{:});
[L,U]=lu(A);
for t = tspan(1)+dt:dt:tspan(2)
    fn1 = feval(f,x,t,varargin{:});
    rhs = An*un+dt*(theta*fn1+(1-theta)*fn);
    temp = feval(g,[xspan(1),xspan(2)],t,varargin{:});
    rhs([1,N]) = temp;
    u = L\rhs;
    u = U\u;
    fn = fn1;
    un = u;
end
return