% SIMPSON
%
% SIMPSON(F,A,B,N)
%  
% devuelve el valor aproximado de la integral de F entre A y B 
% con la regla compuesta de Simpson con N+1 puntos
% 
% F debe estar vectorizada
%
% Extraido del texto:
%
%   "Matlab en cinco lecciones de numerico"
%    por V. Dominguez y M.L. Rapun.
%
% Mas informacion en
%
%   http://www.unavarra.es/personal/victor_dominguez


function s=simpson(f,a,b,n)
n=2*ceil(n/2);                          % n es ahora par
h=(b-a)/n;                              % calculamos h
x=linspace(a,b,n+1);                    % construimos malla
y=feval(f,x);                           % evaluamos f en la malla
s=h/3*(y(1)+...
    4*sum(y(2:2:n))+2*sum(y(3:2:n))+...
    y(n+1));                            % aplicamos la regla
                       
return