Given points , derive an approximation for the integral 
Suppose we are given points $(x_i, f(x_i))$ and $(x_i, f'(x_i)$) for
$i=1,2$. Use Newton Divided differences to design a cubic polynomial $P(x)$
so that
$$ \int\limits_{x_1}^{x_2} P(x) dx = (x_2 - x_1) \frac{ f(x_1) +
 f(x_2) }{2} + \frac{ (x_2-x_1)^2 }{12 } (f'(x_1) - f'(x_2) ) $$

Attempt
We want $P(x_i) = f(x_i)$ and $P'(x_i) = f'(x_i)$. First of all, notice that if we put
$$ H(x) = f[x_1,x_2](x-x_2) + f(x_2) $$
Then $H$ interpolates the first two points. I want to incorporate perhaps adding the term
$$ f'[x_1,x_2](x-x_1) $$
but we need to modify so as to interpolate the other points. Am I on the right track?
 A: Judging by your question, you may have only seen divided differences and the interpolation formula for distinct $x_i$.  The idea is to change the fitting of derivative $f'(x_i)$ to a point $(x_i+\varepsilon, f(x_i)+\varepsilon f'(x_i))$ infinitesimally close to $x_i$ (similarly higher derivatives if needed).
The usual divided difference table for four points looks like
$$
\begin{array}{c|c|c|c|c}
x_0 & f(x_0) &&&\\
&&f[x_0,x_1]&&\\
x_1 & f(x_1) &&f[x_0,x_1,x_2]&\\
&&f[x_1,x_2]&&f[x_0,x_1,x_2,x_3]\\
x_2 & f(x_2) &&f[x_1,x_2,x_3]&\\
&&f[x_2,x_3]&&\\
x_3 & f(x_3) &&&
\end{array}
$$
In our case, we put $x_0$ infinitesimally close to $x_1$ and $x_3$ infinitesimally close to $x_2$.  The divided difference table therefore starts with
$$\require{color}
\begin{array}{c|c|c|c|c}
x_1 & f(x_1) &&&\\
&&f[x_1,x_1]={\color{red}f'(x_1)}&&\\
x_1 & f(x_1) &&\dots&\\
&&\dots&&\dots\\
x_2 & f(x_2) &&\dots&\\
&&f[x_2,x_2]={\color{red}f'(x_2)}&&\\
x_2 & f(x_2) &&&
\end{array}
$$
the $\dots$ you build as usual, e.g.,
\begin{align*}
f[x_1,x_2]&=\frac{f(x_2)-f(x_1)}{x_2-x_1},\\
f[x_1,x_1,x_2]&=\frac{f[x_1,x_2]-f[x_1,x_1]}{x_2-x_1}.
\end{align*}
Now taking limit $x_0\to x_1$, $x_3\to x_2$ in Newton's divided differences interpolation formula
$$
P(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+f[x_0,x_1,x_2,x_3](x-x_0)(x-x_1)(x-x_2)
$$
gives
$$
P(x)=f(x_1)+f[x_1,x_1](x-x_1)+f[x_1,x_1,x_2](x-x_1)^2+f[x_1,x_1,x_2,x_2](x-x_1)^2(x-x_2)
$$
which you can integrate.
