
%-------------------------------------------------------------------------------
% QUESTO CODICE SEGUE [ATKINSON, AN INTRODUCTION TO NUMERICAL ANALYSIS, p.414].
% NECESSITA DELLA FUNCTION: "g".
%-------------------------------------------------------------------------------

demoexample=3;   % SETTA EQUAZIONI E TEMPI DIVERSI. VALORI: 1,2.

theta=0.5;         % [theta=0]   EULERO ESPLICITO. 
                 % [theta=0.5] CRANK-NICOLSON. 
                 % [theta=1]   EULERO IMPLICITO.

mvect=[4 8 16 32 64];  % VETTORE DI "m".
mvect=32;
time_step_factor_vect=ones(size(mvect));                  
                 
switch demoexample
case 1
    
    tfin=5;          % TEMPO FINALE.
    
    G=inline('(-0.1+pi^2)*(exp((-0.1)*t).*sin(pi*x))','x','t');
    d0=inline('zeros(size(t))','t');
    d1=inline('exp((-0.1)*t).*sin(pi)','t');
    f=inline('sin(pi*x)','x');   
    
    solution=inline('exp((-0.1)*t).*sin(pi*x)','x','t');
    
    case 2
        
    tfin=0.1;        % TEMPO FINALE.
    
    G=inline('zeros(size(x))','x','t');
    d0=inline('zeros(size(t))','t');
    d1=inline('zeros(size(t))','t');
    f=inline('sin(pi*x)','x');   
    
    solution=inline('exp((-pi^2)*t).*sin(pi*x)','x','t');
    
    case 3
        
    tfin=1;        % TEMPO FINALE.
    
    G=inline('-sin(x).*sin(t)+sin(x).*cos(t)','x','t');
    d0=inline('zeros(size(t))','t');
    d1=inline('sin(1)*cos(t)','t');
    f=inline('sin(x)','x');   
    
    solution=inline('cos(t).*sin(x)','x','t');
    
end

err_hist_prev_m=[];

fprintf('\n \t [THETA]: %3.3f [TFIN]: %3.3f',theta,tfin);

for mindex=1:length(mvect)
    
    
    m=mvect(mindex); % SELEZIONE m TRA UNA LISTA.
    time_step_factor=time_step_factor_vect(mindex); % TIME FACTOR.
    err_hist=[];
    
    [V_hist,x_mid,t_mid,hx,ht]=cnheat(theta,tfin,m,f,d0,d1,G,time_step_factor);
    
    [X, Y]=meshgrid(x_mid,t_mid);  % VALUTAZIONE ERRORI.  
    U=feval(solution,X,Y);         % SOL. ESATTA NELLA GRIGLIA.
    
    err=norm( U(size(U,1),:)-V_hist(size(V_hist,1),:), inf); % ERRORE COMPIUTO PER "m" FISSATO.
    err_hist=[err_hist; err];
    
    fprintf('\n \t [m]: %3.0f [ERROR]: %2.2e [hx]: %2.2e [ht]: %2.2e', m, err,hx,ht); 
    
    if length(err_hist_prev_m) > 0
        fprintf(' [RATIO]: %2.2f', err_hist_prev_m(size(err_hist_prev_m,1))/err );
    end
    
    err_hist_prev_m=err_hist;    
    
end

mesh(X,Y,U)