% Input: T, matrice m x m triangolare superiore quadrata
%        f, matrice m x k dei termini noti ridotti  
% Output: f, matrice m x k, soluzione  
function f = sostituzione( T , f )
  m = size(T, 1) ;
  % Ciclo sugli m valori incogniti da determinare, all'indietro
  for i = m : -1 : 1
    jx  = (i+1) : m   ; % indici dei termini incogniti già risolti
    v_T = T( i , jx ) ; % corrispondenti termini di T per l'incognita i
    v_f = f( jx , : ) ; % termini incogniti già risolti
    f( i , : ) = ( f( i , : ) - v_T * v_f ) / T( i , i ) ;
  end
end