function s=simpson(f,a,b,N)    

% Calculo de integrales mediante el metodo de Simpson.    
%
%  Entradas:  f = funcion que deseamos integrar 
%             introducida como una cadena de caracteres (entre comas '...').
%         
%    a, b = extremos del intervalo de integracion.    
%       N = numero PAR de subintervalos de la particion de [a,b].
%  
% Salidas:     s =  valor aproximado de la integral.     
%          error =  estimacion del error cometido.
%
%  Comandos de Matlab utilizados (para display inicial y final): 
%
%  disp('...') = display de una cadena de caracteres
%  mat2str(x)  = convierte la "matriz" x en una cadena de caracteres
%
      
disp(' ');
disp('Calculo de la integral mediante el Metodo de Simpson');  
disp(' ');  
disp(['Funcion que se integra: F(x) = ' f ]); 
disp(['Intervalo de integracion: [' mat2str(a) ',' mat2str(b) ']' ]);

       
x=a;   
fa=eval(f);   
x=b;  
fb=eval(f);  
h=(b-a)/N;  
s1=0;   

for i=1:(N/2)
x=a+(2*i-1)*h;
s1=s1+eval(f);   
end   
s2=0;    
for i=1:(N/2-1);
x=a+2*i*h;
s2=s2+eval(f);  
 end  

 s=(h/3)*(fa+4*s1+2*s2+fb);     

disp(' ');
disp(['La solucion aproximada es:' mat2str(s)]); 
disp(' ');

% Estimacion del error 

M= input('Introduce una cota M=sup|f创创|: ' );
error=M*(b-a)*h^4/180;

disp(' ');
disp(['La cota de error es:' mat2str(error)]);  
disp(' ');
   