METODO DEI TRAPEZI COMPOSTO

Ciao a tutti, non capisco cosa non vada nel codice che ho scritto sul metodo dei trapezi composto:
In output S è la somma del metodo dei trapezi, err è la differenza tra le somme di due iterate successive, ierr è un indicatore di errore che vale -1 se è stato superato il numero massimo di iterate consentite nfmax, h è l'ampiezza di ciascun sottointervallo. In input fun è la funzione, a e b sono gli estremi dell'intervallo, tol è la tolleranza.

 function [S, err,ierr,h]=mytrap(fun, a, b,tol,nfmax)
ierr= 0;
err= 3*tol+1; 
k= 1; 
fa= fun(a); 
fb= fun(b);
fab= fa + fb; 
Stemp= 0; % azzero la somma: devo azzerarla una sola volta 
h= b-a; % l'ampiezza dell'intervallo con cui faccio la prima iterata 
AIOLD= fab*h/2; %alla prima iterata ci sono soltanto gli estremi dell'intervallo
while (err > 3*tol && k <= nfmax)
    m= 2^k 
    h= (b-a)/2^k
    
    for i= 1: m % la somma la faccio dal primo al penultimo punto, ho infatti m+1 punti
         Stemp= Stemp + fun(a+ (2*i -1)*h); % alla prima iterata m=2 e i nodi sono 3, quindi i nodi sembrano essere m+1 
    end
    
    S= (2*Stemp + fab)*h/2;
    err = abs(S -AIOLD);
    k= k+1; 
    AIOLD= S; 
end
err= err/3; 
disp ('errore')
disp(err) 
if k> nfmax % indicatore di errore per il numero di iterate 
    disp('superato numero di iterate')
    ierr = -1;
end
  disp('ampiezza intervallo')  

Qualcuno può aiutarmi? il codice non dà errore ma non mi restituisce la somma giusta.