% ---------------------------------------------
% Integral using trapezoidal and simpson rules
% (see also what happens to simpson result if 
%  there are even number of points)  
% ---------------------------------------------
clear all;close all;
format long
exact = 1.6583475942188740493;
f = @(x) sin(x)./x; 
Nbeg = 11;Nend=1000;Nstep=10;

trap=[];
for N=Nbeg:Nstep:Nend
    x   = linspace(0,10,N);
    h   = x(2)-x(1);
    res = 1+f(x(N));   % set manually sin(0)/0 = 1
    w = 2;
    for i=2:N-1
        res = res + w*f(x(i));
    end
    res  = res*h*0.5;
    trap = [trap [h res]'];
end


simp=[];
for N=Nbeg:Nstep:Nend
    x   = linspace(0,10,N);
    h   = x(2)-x(1);
    res = 1+f(x(N));
    w   = 2;
    for i=2:N-1
        w   = 6-w ;
        res = res + w*f(x(i));
    end
    res  = res*h/3;
    simp = [simp [h res]'];
end

% output 
fprintf('%5s %20s %20s \n','N','trap. error','simpson error');
i = 1;
for N=Nbeg:Nstep:Nend
    fprintf('%5i %20.15f %20.15f\n',N,trap(2,i)-exact,simp(2,i)-exact);
    i = i + 1;
end

plot(trap(1,:),abs(trap(2,:)-exact),'b-',simp(1,:),abs(simp(2,:)-exact),'g-');
legend('trap','simp');
xlabel('stepsize h');ylabel('|error|')
    