% ---------------------------------
% Integral using trapezoidal rule
% ---------------------------------
clear all;close all;
format long
exact = 1.6583475942188740493;
f = @(x) sin(x)./x; 
fprintf('%5s %20s %20s %20s\n','N','h','res','res-exact');
data=[];
for N=11:10:1000
    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;
    fprintf('%5i %20.15f %20.15f %20.15f\n',N,h,res,res-exact);
    data = [data [h res-exact]'];
end
plot(data(1,:),abs(data(2,:)),'b-')
xlabel('stepsize h');ylabel('|error|')


