(* Verify Simpson' s rule *)
(* Use build in Interpolating polynomial. *)

Clear[a, b];
data = {{a, f0}, {(b + a)/2, f1}, {b, 
            f2}};
f = InterpolatingPolynomial[data, x];
intf = Simplify[Integrate[f, {x, a, b}]]

(* Verify Simpson 3/8 rule *)

!! This cell is written for reminding the formula!!
(∫_x1^x4) 
          f (x) \[DifferentialD]x = (3/8) (f[x1] + 3  f[x2] + 
              3  f[x3] + f[x4]) Δx  +  
        O (Δ(x^5) 
              (f^((4)))) 
F[x] = (∫_x1^x) f[x] \[DifferentialD]x
D[F[x], x] = f[x]

Clear[a, b];
Δx = (b - a)/3;
data = {{a, f0}, {a + Δx, 
            f1}, {a + 2 Δx, f2}, {b, 
            f3}};
f = InterpolatingPolynomial[data, x];
intf = Simplify[Integrate[f, {x, a, b}]]

(* Use Taylor expansion. *)

Clear[f, F, fpoly, Fpoly, x]
fpoly = Series[f[x], {x, x1, 5}]
Fpoly = Series[F[x], {x, x1, 6}]

substitution={F'[x1]->f[x1],F''[x1]->f'[x1],Derivative[3][F][x1]-> f''[x1],Derivative[4][F][x1]->Derivative[3][f][x1],Derivative[5][F][x1]->Derivative[4][f][x1]};

Clear[fpoly, Fpoly]
fpoly = Normal[Series[f[x], {x, x1, 4}]]
Fpoly = Normal[Series[F[x], {x, x1, 5}]] /. substitution

x2 = x1  +  Δx;
x3 = x1  +  2 Δx;
x4 = x1  +  3 Δx;
f1 = f[x1];
f2 = fpoly /. x -> x2;
f3 = fpoly /. x -> x3;
f4 = fpoly /. x -> x4;
F1 = F[x1];
F4 = Fpoly /. x -> x4;
Simplify[
      F4 - F1 - (3/8) (f1 + 3 f2 + 3  f3 + 
              f4) Δx]