function [ret, x] = simpson (fun, n, a, b)
%
% [ret, x] = simpson (fun, n, a, b)
%
% Formula di Simpson a n+1 punti (n pari)
h = (b - a) / n;
w = [1, repmat([4, 2], 1, n / 2 - 1), 4, 1];
x = linspace (a, b, n + 1)';
ret = h / 3 * (w * fun (x));
%!demo
%! fun = @(x) exp (x);
%! a = 0;
%! b = 1;
%! nrange = 10:10:50;
%! counter = 0;
%! for n = nrange
%!   counter = counter + 1;
%!   ret = simpson (fun, n, a, b);
%!   err(counter) = abs (exp(1) - 1 - ret);
%! end
%! loglog (nrange, err, '*', nrange, nrange .^ (-4))
