function I = simpson(fun,a,b,m)

% I = simpson(fun,a,b,m)
%
% Applica il metodo di Simpson 
% per calcolare una approssimazione dell'integrale 
% tra a e b di fun suddividendo [a,b] in 
% m sottointervalli di uguale ampiezza.
% Si assume (b>=a).

h = (b-a)/m; % ampiezza dei sottointervalli

xestremi = a + h*[0:m];
xinterni = a + h/2 +[0:m-1]*h;

festremi = feval(fun,xestremi);
finterni = feval(fun,xinterni);
I = (h/6)*(festremi(1) + 4*sum(finterni) + ...
           2*sum(festremi(2:end-1)) + festremi(end));
       
return

