This problem reduces to linear algebra. Start with
$$\int_0^1ug(u)du=ag(0)+bg^{\prime}(0)+cg(1)+dg^{\prime}(1)$$
for $g(u)=\in\{1,u,u^2,u^3\}$ and you get
$$\begin{align}\frac12&=a+c\tag{1}\\
\frac13&=b+c+d\tag{2}\\
\frac14&=c+2d\tag{3}\\
\frac15&=c+3d\tag{4}\end{align}$$
Subtracting eq $(3)$ from eq $(4)$, we find that $d=-\frac1{20}$. Substituting into eq$(3)$, we have $c=\frac7{20}$. Then substituting into eq $(2)$ and $(1)$ we show that $b=\frac1{30}$ and $a=\frac3{20}$. Now we transform the original integral via $\frac{x-x_0}{x_1-x_0}=u$; $x=x_0+(x_1-x_0)u$
$$\begin{align}\int_{x_0}^{x_1}(x-x_0)f(x)dx&=\int_0^1(x_1-x_0)uf(x_0+(x_1-x_0)u)(x_1-x_0)du\\
&=(x_1-x_0)^2\int_0^1uf(x_0+(x_1-x_0)u)du\\
&=(x_1-x_0)^2\left[\frac3{20}f(x_0)+\frac1{30}(x_1-x_0)f^{\prime}(x_0)+\frac7{20}f(x_1)-\frac1{20}(x_1-x_0)f^{\prime}(x_1)\right]\end{align}$$
Comparing with the original requirements,
$A=\frac3{20}$, $B=\frac7{20}$, $C=\frac1{30}$, and $D=-\frac1{20}$. Checking, we find that the solution is indeed exact for $f(x)$ a polynomial of degree at most $3$.
Forgot the remainder: since the formula is interpolatory and there are no interior nodes, we know the remainder has to have the form $R=Ch^6f^{(4)}(\xi)$ for some $\xi\in(x_0,x_1)$. We can find $C$ by attempting $f(x)=(x-x_0)^4$ with our formula, and we find
$$\frac16h^6=\frac7{20}h^6-\frac1{20}\cdot4h^6+R=\frac3{20}h^6+Ch^6\cdot24$$
So $C=\frac1{1440}$ and $R=\frac1{1440}h^6f^{(4)}(\xi)$.
EDIT: The form for the remainder can be derived in a fashion similar to the way the remainder term for the Taylor series is derived. We know that the quadrature formula is equivalent to finding a cubic polynomial $p(u)$ that matches the function $g(u)$ and its first derivative for $u\in\{0,1\}$. But now, because we have nothing better to do, we want to find a quartic polynomial $q(u)$ that does the same and also $q(c)=g(u)$ for a prescribed $c\in(0,1)$. It's pretty easy to do, given $p(u)$:
$$g(u)=p(u)+\frac{u^2(u-1)^2}{c^2(c-1)^2}\left(g(c)-p(c)\right)$$
Just check that it matches all $5$ conditions! Now we want to analyze the deviation of $q(u)$ from $g(u)$. We call that deviation $e(u)=g(u)-q(u)$. Becaue of $3$ of the conditions imposed on $q(u)$, we know that $e(u)=0$ for $u\in\{0,c,1\}$. Because $u=0$ and $u=1$ were double roots of $e(u)$, we know that they are also zeros of $e^{\prime}(u)$. By Rolle's theorem, we know also that there is some $a\in(0,c)$ such that $e^{\prime}(a)=0$ and some $b\in(c,1)$ such that $e^{\prime}(b)=0$. Thus $e^{\prime}(u)$ has $4$ zeros in $[0,1]$. Applying Rolle's theorem $3$ more times we can show that there is some $\eta\in(0,1)$ such that $e^{(4)}(\eta)=0$. This expands to
$$e^{(4)}(\eta)=g^{(4)}(\eta)-q^{(4)}(\eta)=g^{(4)}(\eta)-\frac{24}{c^2(c-1)^2}\left(q(c)-p(c)\right)=0$$
We can turn that around to say that for any $u\in(0,1)$,
$$g(u)=p(u)+\frac{u^2(u-1)^2}{24}g^{(4)}(\eta)$$
For some $\eta\in(0,1)$. Now we are ready to get that remainder!
$$\int_0^1ug(u)du=\int_0^1up(u)du+\frac1{24}\int_0^1u^3(u-1)^2y^{(4)}(\eta(u))du$$
Where we don't know the exact value of $\eta(u)$, just that it lies in $(0,1)$. Comparing with our definition of the remainder we can see that
$$R=\frac1{24}\int_0^1u^3(u-1)^2y^{(4)}(\eta(u))du$$
Since $u^3(u-1)^2\ge0$ for $u\in(0,1)$, the biggest $u^3(u-1)^2y^{(4)}(\eta(u))$ could be in $(0,1)$ is $u^3(u-1)^2y^{(4)}(\eta(u))\le u^3(u-1)^2\max\left(y^{(4)}(\eta(u))\right)$. Thus
$$R\le\frac1{24}\max\left(y^{(4)}(\eta)\right)\int_0^1u^3(u-1)^2du$$
Similar arguments apply at the lower bound so that
$$\frac1{24}\min\left(y^{(4)}(\eta)\right)\int_0^1u^3(u-1)^2du\le R\le\frac1{24}\max\left(y^{(4)}(\eta)\right)\int_0^1u^3(u-1)^2du$$
Given continuous fourth derivative, it follows by the intermediate value theorem that
$$R=\frac1{24}y^{(4)}(\xi)\int_0^1u^3(u-1)^2du$$
For some $\xi\in(0,1)$. We can knock down that integral to show that $R=\frac1{1440}y^{(4)}(\xi)$. The only thing that could have stopped our theorem from working for smooth functions was the requirement that $u^3(u-1)^2$ have the same sign for $u\in(0,1)$. Thus, if the weight function (here $(x-x_0)\rightarrow u$) doesn't change sign in the interval of integration and any nodes lying within the interval of integration have their first derivatives either incorporated in the integration formula or are ignorable (as for Gaussian integration formulas) the ansatz is going to be valid. I posted an example where these conditions are not met and the ansatz is invalid. But that's not a problem for us. The factor of $h^6$ in the ultimate formula for the remainder has contributions of $h^2$ from the transformation from $[x_0,x_1]$ to $[0,1]$ and of $h^4$ from the $4$ derivatives taken.
\begin{align}1&=a+c\tag{1}\ \frac12&=b+c+d\tag{2}\ \frac13&=c+2d\tag{3}\ \frac14&=c+3d\tag{4}\end{align}
right?
– Relure May 16 '16 at 21:43