function area=simpson(fun,a,b,m)
%Aproxima por la regla de los trapecios compuesta el valor de la integral 
%de una función fun(x) en un intervalo de extremos a y b tomando m+1 puntos equiespaciados. 
%Variables de entrada:
%     fun(x): funcion que se quiere integrar y que debe 
%     introducirse con notación simbolica (eg. 'g').
%     a: extremo izquierdo del intervalo
%     b: extremo derecho del intervalo 
%     2m: número de puntos memos uno o número de subintervalos.    
% Variables de salida:
%     area: integral aproximada

h=(b-a)/(2*m);
x=a:h:b;
sum_pares=0;
sum_impares=0;
for i=1:m
   sum_pares=sum_pares+feval(fun,x(2*i));
end
for i=1:m-1
   sum_impares=sum_impares+feval(fun,x(2*i+1)); 
end
area=(h/3)*(feval(fun,a)+feval(fun,b)+4*sum_pares + 2*sum_impares);