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

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

% first compute area with n-1 intervals
[area1,h,err] = trap(a,b,n);

% now compute area with 2*(n-1) intervals which means 2*(n-1)+1 meshpoints:
[area2,h,err] = trap(a,b,2*(n-1)+1);

area = (4*area2-area1)/3;
err = (f_int(b)-f_int(a))-area;