Here are some hints.
The approximations that you want are called open Newton-Cotes formulas.
They're constructed in the same way as the closed NC formulas (which are available in Maple directly in Student:-Calculus1:-ApproximateInt), except the endpoints a and b are not used.
So, construct a list of data points [[a+h, f(a+h)], ..., [b-h, f(b-h)]].
Create the interpolating polynomial with CurveFitting:-PolynomialInterpolation.
Integrate the polynomial over the segment [a,b], simplify (alternatively, think about how to make everything simpler from the start by using integers instead of a+h, a+2h, ...), and you're done.
To estimate the error term, apply your code to f on the interval [0, h], subtract int(f(x),x=0..h) and expand the result into a series.