% numerical integration of (1 + 4x^2)^-1 on [-1, 1]

f=@(x) (1+4.*x.^2).^(-1);
errors=zeros(40);
sum=0;

%trapezoid rule
for n=1:40;
   xnodes = linspace(-1,1,n+1); %sets up x nodes
   ynodes=f(xnodes);
   sum=0;
   for i=1:n;
     sum=sum+(ynodes(i)+ynodes(i+1))/n;
   end
   errors(n)=abs(atan(2)-sum); 
end

figure;
loglog(1:40,errors);
xlabel('n');
ylabel('error');
print -depsc 3a.eps;


errors=zeros(40);
%Newton-Cotes quad.
for n=1:40;
   x = linspace(-1,1,n+1); % n+1 equally-spaced x values in [-1,1]
   A = vander(x);           
   A = A(:,end:-1:1); % reverse order of cols
   A=A'; %sets up matrix correctly
   b=zeros(n+1,1);
   for i=1:n+1
      b(i)=1/(i)*((1)^(i)-(-1)^(i)); %sets up b matrix;
   end
   w=A\b; %solves for weights
   
   errors(n)=abs(atan(2)-f(x)*w);
   
end

figure;
loglog(1:40,errors);
xlabel('n');
ylabel('error');
print -depsc 3b.eps;