function [s weights nodes]=numint(func,a,b,n,method)

% Numerically integrate a function, called via a MATLAB mfile func,
% from a to b, using n nodes, using
%
%   sum(weight_i*f(node_i))
%
% where the function evaluation f(node_i) is calculed by calling
% feval(func(nodes))
%
% Method to calculate the weights and nodes:
%   1. Riemmann sums
%   2. Trapezoid Rule
%   3. Simpsons Rule
%   4. Gaussian Quadrature
%
% References for the forst three can be found in any decent calculus
% textbook.   For Gauss Quadrature via the Lanczos matrix, see Trefethen
% and Bau, Numerical Linear Algebra, 1997, or Bai, Fahey and Golub 
% "Some Large Scale matrix computation problems"  in the Journal of Comp. 
% and Appl. Math, 1996.



weights=ones(1,n);
tic
if method==1  % Riemann sums
    weights=(b-a)/n*weights;
    nodes=linspace(a+(b-a)/n,b,n);
    titlestr='Riemann sums';
   
elseif method==2 % Trapezoid Rule
   weights=(b-a)/(n-1)*weights;
   weights([1 end])=weights([1 end])/2;
   nodes=linspace(a,b,n); 
   titlestr='Trapezoid Rule';
   
elseif method==3  % Simpsons Rule
   if ~mod(n,2)
       error('For Simpsons Rule, n, the number of nodes, must be odd')
   end
   weights=(b-a)/(3*(n-1))*weights;
   weights(2:2:n)=4*weights(2:2:n);
   weights(3:2:n-1)=2*weights(3:2:n-1);
   nodes=linspace(a,b,n);
   titlestr='Simpson quadratics';
   
elseif method==4  % Gaussian Quadrature
   dt=.5*1./sqrt(1-1./(2*(1:(n-1))).^2);
   T=spdiags([[dt 0]' [0 dt]'],[-1 1],n,n); 
   [W nodes]=eig(full(T));
   %nodes=diag(weights)';
   %weights=2*W(1,:)'.^2;
   
   % Need to change of variables for a and b other than -1 and 1
   nodes=(b-a)/2*diag(nodes)' + (a+b)/2;
   weights=(b-a)*W(1,:).^2;
   titlestr='Gaussian Quadrature via Lanczos Tridiagonal';
end

    


s=feval(func,nodes)*weights';
toc

%keyboard

% Plot the function over the interval [a, b]
figure(1)
clf
xvals = linspace(a,b,100);
plot(xvals,feval(func,xvals))
hold on
plot(nodes,zeros(length(nodes),1),'r.')
title(sprintf('Numerical Integration using %s',titlestr))
legend(sprintf('%s',func),'nodes')