function [V_hist,x_mid,t_mid,hx,ht]=cnheat(theta,tfin,m,f,d0,d1,G,sigma)
% INPUT.
% theta   : PARAMETRO DEL THETA METODO (TRA 0 ED 1, POSSIBILMENTE 0,1 INCLUSI).
% tfin    :  ISTANTE FINALE.
% m       : DETERMINA IL PASSO SPAZIALE CON PASSO h_x=1/m.
% f       : u(x,0)=f(x)
% d0, d1  : u(0,t)=d0(t), u(1,t)=d1(t).
% G       : (D_t)u=(D^2_x)u+G.
% sigma: IL PASSO TEMPORALE E' sigma*(hx^2)/2.
% OUTPUT.
% V_hist  : SOLUZIONE ISTANTE PER ISTANTE. RIGHE INIZIALI -> BASSO t.
% x_mid   : PUNTI INTERNI NELLA VARIABILE x.
% t_mid   : PUNTI INTERNI NELLA VARIABILE t.
% hx      : PASSO SPAZIALE.
% ht      : PASSO TEMPORALE.

hx=1/m; % DELTA=hx
ht=sigma*(hx^2)/2;    % STEP TEMPORALE.
x=(0:hx:1)'; x_mid=x(2:end-1,:); % PUNTI INTERNI.
u0=feval(f,x_mid);               % VETTORE SOLUZIONE AL TEMPO "0".
t=(0:ht:tfin)';  t_mid=t(2:end,1);   % TEMPI DA ANALIZZARE.
lambda_matrix=(1/hx^2)*(diag(-2*ones(m-1,1))+diag(ones(m-2,1),1)+...
    diag(ones(m-2,1),-1)); % COSTRUZIONE MATRICE LAMBDA.

if theta < 1, gt_prev=g(0,x_mid,hx,d0,d1,G); end 
V_old=u0; V_hist=[V_old']; err_hist=[]; % INIZIALIZZAZIONI.

for index=2:length(t) % STUDIO SOLUZIONE IN TEMPI FINO A tfin.
    % STUDIANDO u'(t_curr)=LAMBDA*u(t_curr)+g(t_curr),
    % DOVE g(t)=(1/hx^2)*(d0(t),0,...,0,d1(t))+[G(x1,t),...,G(x_(m-1),t)].
    t_curr=t(index); gt_curr=g(t_curr,x_mid,hx,d0,d1,G);
    switch theta
        case 0 % EULERO ESPLICITO.
            V_new=V_old+ht*(lambda_matrix*V_old+gt_prev);
            gt_prev=gt_curr;
        case 1 % EULERO IMPLICITO.
            A=eye(size(lambda_matrix))-ht*lambda_matrix;
            b=V_old+ht*gt_curr; V_new=A\b;
        otherwise % THETA METODO.
            A=eye(size(lambda_matrix))-(ht*theta)*lambda_matrix;
            b=V_old+(ht*(1-theta))*lambda_matrix*V_old+...
                (ht*(1-theta))*gt_prev+(ht*theta)*gt_curr;
            V_new=A\b; gt_prev=gt_curr;
    end
    V_hist=[V_hist; V_new']; V_old=V_new;
end
