Instead of using the above first-order polynomials approximations for $f$, we attempt at using a second-order polynomials. In this case we need three points in order to define a second-order polynomial approximation $$\begin{equation*} f(x) \approx P_2(x)=a_0+a_1x+a_2x^2. \end{equation*}$$ Using again Lagrange's interpolation formula we have $$\begin{equation*} P_2(x)=\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}y_2+ \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}y_1+ \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}y_0. \end{equation*}$$ Inserting this formula in the integral of Eq. (2) we obtain $$\begin{equation*} \int_{-h}^{+h}f(x)dx=\frac{h}{3}\left(f_h + 4f_0 + f_{-h}\right)+O(h^5), \end{equation*}$$ which is Simpson's rule.