function [q, tt]=edromberg(a,b,n);
% EDROMBERG by Ed Bueler 11/6/02
% Approximates the integral from a to b of f(x)dx by using
% the trapezoidal rule with Richardson extrapolation.
% Trap. rule is used with increasing numbers of subdivisions,
% up to 2^n. 
% *Inefficiently* computes the trapezoid approximations.
% *Inefficiently* uses polyfit,polyval to do the extrapolation.

for j=1:n
  hh(j)=(b-a)/(2^j);
  tt(j)=trap(a,b,2^j);
end
q=polyval(polyfit(hh.^2,tt,n-1),0);
tt=tt';

function t=trap(a,b,k);
% helper function to EDROMBERG
% performs trapezoid rule with k subintervals

x=linspace(a,b,k+1); % a row vector
h=(b-a)/k;
y=f(x');  % f defined below; y is a column vector
ct=[1 repmat(2,1,k-1) 1]; % coeffs for trap; a row vector
t=.5*h*ct*y;


function y=f(x);
% defines integrand

y=(1+x.^4).^(-1);
