% ---------------------------------
% Integral using trapezoidal rule
% ---------------------------------
clear all;close all;
format long
exact = 1.6583475942188740493;
f = @(x) sin(x)./x; 
disp('h, results and error')
out=[];
for N=21:50:200
    x=linspace(0,10,N);
    h = x(2)-x(1);
    res = 1+f(x(N));
    w = 2;
    for i=2:N-1
        res = res + w*f(x(i));
    end
    res = res*h*0.5;
    out = [out ; [h res res-exact]];
end
disp(out);
    