%NC demo
clc
clear all
close all
%Newton Cotes integration of the cosine function
%Range is 0 to pi/2
%fname = 'humps';
%fname = vectorize(inline('cos(x)-1/sqrt(x+0.01)'));
fname = 'cos';
a=0;
b=pi/2;
z = linspace(a,b,100);
mlim=11;
hold on
fexact=feval(fname,z);
exactinteg=quad(fname,a,b,1.0e-8);
plot(z,fexact)
for m = 2:1:mlim
    c = WNC(m);
    fx = zeros(1,m);
    x = linspace(a,b,m);
    fx = feval(fname,x);
    integ=(b-a)*(fx*c);
    error=abs(integ-exactinteg);
    disp (sprintf(' %2.0f points Int= %20.16f Err=%20.7e',m,integ,error))
    acoef=polyfit(x,fx,m-1);
    fz = polyval(acoef,z);
    figure
    plot(z,fz,z,fexact,'-r',x,fx,'*')
    pause
end
hold off