function tr = trapzd(f,rng,N)

a = rng(1); b = rng(2);
h = (b-a)/N;
x = a + h*(1:(N-1));
tr = h*(0.5*(f(a)+f(b))+sum(f(x)));