Let $F$ be an anti-derivative of $f$, $F'=f$. W.l.o.g. $x_1=0$, set $x=h$ then we are interested in the error expression
$$
g(x)=F(x)-F(-x)-\frac{x}3(f(x)+4f(0)+f(-x)).
$$
This has derivatives
\begin{alignat}{2}
g'(x)&=\frac23(f(x)-2f(0)+f(-x))&&-\frac x3 (f'(x)-f'(-x))
\\
g''(x)&=\frac13(f'(x)-f'(-x))&&-\frac x3(f''(x)+f''(-x))
\\
g'''(x)&=&&-\frac x3(f'''(x)-f'''(-x)).
\end{alignat}
Consequently, by the extended mean value theorem
$$
\frac{g(x)}{x^m}
=\frac{g'(x_1)}{mx_1^{m-1}}
=\frac{g''(x_2)}{m(m-1)x_2^{m-2}}
=\frac {g'''(x_3)}{m(m-1)(m-2)x_3^{m-3}}
$$
with $0<x_3<x_2<x_1<x$.
Using $m=5$ this gives
$$
\frac{g(x)}{x^5}=\frac{g'''(x_3)}{60x_3^2}=-\frac1{90}·\frac{f'''(x_3)-f'''(-x_3)}{2x_3}=-\frac1{90}·f^{(4)}(x_4).
$$
with $|x_4|<x_3<x$.
This results in the error formula
$$
g(x)=-\frac1{90}·f^{(4)}(x_4)·x^5.
$$
or after translating the initial simplifications back,
\begin{align}
\int_{x_0}^{x_2}f(x)\,dx
&=\frac{h}{3}[f(x_0)+4f(x_1)+f(x_2)]-\frac{h^5}{90}f^{(4)}(\xi)
\end{align}
with $\xi\in(x_0,x_2)$
Once the the quadrature formula obtained by cubic interpolation of $f$ based at nodes $(a,f(a))$, $(a+h,f(a+h))$, $(a+2h,f(a+2+2h))$ and $(b,f(b))$ with $h=\frac{b-a}{3}$ and its error are established, notice that this formula can be rewritten as
$$\int^b_af\approx\frac{3h}{8}\Big(f(a)+3f(a+h)+3f(a+2h)+f(b)\Big)$$
with error $e_h=-\frac{(b-a)^5}{6480}f^{(iv)}(c)$ for some $c\in (a, b)$.
Divide the interval $[a,b]$ in $3m$ pieces of same length $h=\frac{b-a}{3m}$. This defines a partition $x_k=a+\frac{b-a}{3m}$.
On the each subinterval of the form $[x_{3k},x_{3(k+1)}]$, $k=0,\ldots ,m-1$ apply Simpson's rules based on cubic interpolation (the fist formula in the OP). That is
$$\begin{align}
\int^b_af=\sum^{m-1}_{k=0}\int^{x_{3(k+1)}}_{x_{3k}}f\tag{1}\label{one}
\end{align}$$
and
$$I_k:=\int^{x_{3(k+1)}}_{x_{3k}}f\approx\frac{3h}{8}\Big(f(x_{3k})+3f(x_{3k+1})+3f(x_{3k+2})+f(x_{3(k+1)})\Big)=A_k$$
Denoting the righthand side of \eqref{one} by $A$ we have that
\begin{align}
\frac{8}{3h}A&=\sum^{m-1}_{k=0}\Big(f(x_{3k})+3f(x_{3k+1})+3f(x_{3k+2})+f(x_{3(k+1)})\Big)\\
&=f(x_0)+\sum^{m-1}_{k=0}3\big(f(x_{3k+1})+f(x_{3k+2})\big) +2\sum^{m-1}_{k=1}f(x_{3k}) + f(x_{3m})
\end{align}
The composite formula follows.
As for the error, each approximation $A_k$ to $I_k$ has error $-\frac{3}{80} h^5 f^{(iv)}(c_k)$ for some $c_k\in (x_{3k},x_{3(k+1)})$. The sum of all errors is
$$E:=-\frac{3}{80}h^5\sum^{m-1}_{j=0}f^{(iv)}(c_k)$$
If we assume that $f\in\mathcal{C}^4([a,b])$, there there is $c\in(a, b)$ such that $f^{(iv)}(c)=\frac1m\sum^{m-1}_{k=0}f^{(iv)}(c_k)$ and so
$$
E=-\frac{1}{80}(b-a)h^4f^{(iv)}(c)
$$
With regards to the error in this quadrature, there are somewhat explicit formulas for the errors resulting from interpolation methods on evenly distributed nodes:
For $n\in\mathbb{N}$, $n\geq2$, the error of approximation to $\int^b_af$ in the quadrature method based on interpolating over the uniform partition $a_k=\frac{k}{n}(b-a)$, $0\leq k\leq n$, is
- If $n$ is even and $f\in C^{(n+2)}[a,b]$,
$$E_n=\frac{K_n}{(n+2)!}f^{(n+2)}(\eta),\quad K_n=\int^b_a x\prod^n_{k=0}(x-x_k)\,dx$$
for some $a<\eta<b$,
- If $n$ is odd and $f\in C^{(n+1)}[a,b]$,
$$E_n=\frac{K_n}{(n+1)!}f^{(n+1)}(\eta),\qquad K_n=\int^b_a\prod^n_{j=1}(x-x_j)\,dx$$
For Simpson's 3/8 rule, $n=3$ and so
$K_3=\int^b_a (x-a)(x-\tfrac{2a+
b}{3})(x-\tfrac{a+2b}{3})(x-b)\,dx$
The simplification, though messy with pencil and paper, yields
$E_3=-\frac{(b-a)^5}{6480}f^{(4)}(\eta)$ for some $a<\eta<b$.
Reference: Krylov, V. I. (translated by Stroud, A.), Approximate calculation of integrals, Dover Publications. 2005 pp 89-91
Best Answer
Since $h=\dfrac {b-a}n$ and $a=0, b=3$:
$$\int_a^b{f(x)dx}\approx S_n = \frac{1}{3}h[f(a)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+...+4f(x_{n-1})+f(b)]\\=\frac 16 ( 3+ 11.8321+ 5.6568+10.3923+4.4721+6.6332+0)\\=6.9977$$