% [h,area,err] = simp(a,b,n)

function [h,area,err] = simp(a,b,n)

h = (b-a)/(n-1);
for i=1:n
  x(i) = a + (i-1)*h;
end

area = 0;
for i=1:(n-1)/2
  k = 2*i;
  area = area + f(x(k-1)) + 4*f(x(k)) + f(x(k+1));
end
area = area*h/3;

area = area - h^4/180*(fppp(b)-fppp(a));


err = area - (f_int(b)-f_int(a));