To recap, the 3-point Simpson rule divides the interval of
integration (x0, x1, x2) into 2 equal panels of width h, and the
3 ordinates at the panel boundaries are combined using simple
arithmetic |S = (h/3) (y0 + 4y1 + y2)| to approximate the actual
integral I. Next, the rule is applied to the two subintervals
(x0, x1) and (x1, x2), and the sum of these results, S next,
i.e., the new best value, is checked to see if it is “close
enough” to S, the most recent old best value. If so, S next is
returned. If not, further subdivision is called for. This is the
so-called “adaptive quadrature”, as the frenzy of subdividing
adapts to those sections of the interval that need it most.
(Note that when subdividing the interval (x0, x1, x2) into (x0,
x0.5, x1) and (x1, x1.5, x2), two new ordinates are needed at
x0.5 and x1.5, but the old ordinates at x0, x1, and x2 can be
reused.) Newton-Cotes formulas with more points give better
results. Simpson’s rule (3-point) has an error term that is
proportial to the 5th power of t