%%%%%%%%%  Newton_Cotes Formula  %%%%%%%%%%
fprintf('\n Newton-Cotes Formula\n');
n=input('Number of data points(2-8)?');
a=input('Lower limit of integration:');
b=input('Upper limit of integration:');
%%%%%%%%%%%%%%%%%%  柯特斯系數  %%%%%%%%%%%%%%%%%%%%%
a_{1}=[1/2,1/2];
a_{2}=[1/6 2/3 1/6];
a_{3}=[1/8 3/8 3/8 1/8];
a_{4}=[7/90 16/45 2/15 16/45 7/90];
a_{5}=[19/288 25/96 25/144 25/144 25/96 19/288];
a_{6}=[41/840 9/35 9/280 34/105 9/280 9/35 41/840];
a_{7}=[751/17280 3577/17280 1323/17280 2989/17280 2989/17280 1323/17280 3577/17280 751/17280];
a_{8}=[989/28350 5888/28350 -928/28350 10496/28350 -4540/28350 10496/28350 -928/28350 5888/28350 989/28350];

n=n-1;
h=(b-a)/n;
x=[a:h:b];
y=f(x);
Result=(b-a)*sum(a_{n}.*y);

%%%%%%%%%%%%%%%  輸出結果  %%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n---------------------------------------------------------\n');
fprintf(' n          x              f(x)              w        \n');
fprintf('---------------------------------------------------------\n');

for i=1:n+1
     fprintf(' %d   %12.5f     %12.5f    %12.5f \n',i-1,x(i),y(i),a_{n}(i));
end

fprintf('---------------------------------------------------------\n');
fprintf('Result    I=%6.5f',Result);
fprintf('\n---------------------------------------------------------\n');