%Algoritmo de Simpson para o calculo de Integrais:
%Calcula o numero de intervalos para que a solucao seja
%identica a solucao numerica do Matlab com uma certa tolernacia
%%SIMPSON Simpson's rule integration with equally spaced points
%
%  [Integral,m]=SIMPSON(f,a,b,n) returns the Simpson's rule approximation to
%  the integral of f(x) over the interval [a,b] using n+1 equally
%  spaced points.  The input variable f is a string containing the
%  name of a function of one variable.  The function f(x) must accept
%  a vector argument and return the vector of values of the function.
%m numero deintervalos para soluç~ao identica ao quad do matlab.
%
%  NOTE: n must be even.

function [Integral,m]=simpson(func,a,b,n)

if rem(n,2)~=0
    beep;
    error('n deve ser par');
end

h=(b-a)/n;

x=a:h:b;
 y=feval(func,x);
%y=sin(x);
for i=1:n+1
    if i==1  | i==n+1
        aux(i)=1;
    else
        if rem(i,2)==0
            aux(i)=4;
        else
            aux(i)=2;
        end
    end
end

Integral=h/3*(sum(y.*aux));
I=quadl(func,a,b);
    m=n;
    while abs((Integral-I))>0.00001
        m=m+2;
        
        h2=(b-a)/m;

x=a:h2:b;
y=feval(func,x);
% y=sin(x);
for j=1:m+1
    if j==1  | j==m+1
        aux(j)=1;
    else
        if rem(j,2)==0
            aux(j)=4;
        else
            aux(j)=2;
        end
    end
 end   
Integral=h2/3*(sum(y.*aux));


end