function ShowNCIdea(m)
% Illustrates the idea behind the Newton-Cotes approach
% to quadrature. The m-point rule is applied to the problem
% of integrating f(x) = exp(-x).*sin(7*x) across [0,3].

close all
a = 0;
b = 3;
% Generate f-values for display
x = linspace(a,b)';
y = f(x);
% Generate interpolation points...
x0 = linspace(a,b,m)';
y0 = f(x0);
% The interpolating polynomial...
p = polyfit(x0,y0,length(x0)-1);
pvals = polyval(p,x);

% Plot the function.
subplot(2,2,1)
plot(x,y,[a b],[0 0],':k')
title('f(x) = exp(-x/2)sin(2*x)','fontsize',14);
axis([a b -1 1])

% Plot the interpolant.
subplot(2,2,2)
plot(x,y,x,pvals,'r',x0,y0,'*r',[a b],[0 0],':k')
title('The Interpolant','fontsize',14);
xlabel(sprintf('m = %1d',m),'fontsize',14)
axis([a b -1 1])

% Display the integral of the function.
subplot(2,2,3)
fill([x;b;a],[y;0;0],[.5 .5 .5])
title('Integral of f(x)','fontsize',14);
axis([a b -1 1])

% Display the integral of the interpolant
subplot(2,2,4)
plot(x,pvals)
fill([x;b;a],[pvals;0;0],'r')
title('Integral of Interpolant','fontsize',14);
axis([a b -1 1])

integralQUAD = quad(@f,a,b,10^(-12));
integralQNC = QNC(@f,a,b,m);
error = abs(integralQUAD-integralQNC);
xlabel(sprintf('Error = %8.3e',error),'fontsize',14)

function y = f(x)
y = exp(-x/2).*sin(2*x);