The formula that you give for the absolute value of the error is exact, sort of. By that I mean that for every $x$, there is a $\xi(x)$ in our interval such that the error when you use the interpolating quadratic at $x$ is exactly the one obtained by putting $\xi=\xi(x)$ in the error formula.
In most cases, this kind of Mean Value Theorem information leads to estimates, but not to exact values, since we don't know $\xi(x)$.
However, here we are evaluating the third derivative of our cubic at $\xi(x)$. The third derivative is constant! So it doesn't matter what $\xi(x)$ is, that part of the error estimate does not change. The rest is an explicit polynomial in $x$. It follows that the error estimate and the actual error are equal, as you noticed computationally. Both are polynomials in $x$. Given any two polynomials, equality at all values of $x$ (or even at infinitely many values of $x$, or, for cubics, at $4$ values of $x$) means that the coefficients match.
Precisely the same thing happens when you use the same process to approximate a polynomial of degree $n+1$ by a Newton interpolating polynomial of degree $n$. If you look up the error estimate,
you will see that it is a certain polynomial in $x$ times a term $f^{(n+1)}(\xi(x))$. (I like that version better than the absolute value version, since we also get information about the sign of the error.)
For a polynomial $f$ of degree $n+1$, the $(n+1)$-th derivative is constant, so
the term $f^{(n+1)}(\xi(x))$ is a constant independent of $x$. Thus the error estimate and the error are identically equal. Since they are both polynomials, that means their coefficients are equal.
To know why Gauss quadrature works, you should look at the proof.
As for how to do it, you need to do Gram-Schmidt on the standard polynomial basis to get a degree three orthogonal polynomial. The roots of this polynomial will constitute the nodes. Then, test your rule on the polynomials $1,x$, and $x^2$ to get a linear system of equations to find the weights. The proof explains why this is what you want to do, and this is the algorithm that you always follow for such a problem.
EDIT: How it's done-
First, we apply Gram-Schmidt to the standard basis $v_1=1,v_2=x,v_3=x^2,v_4=x^3$ with respect to the inner product $\langle f,g\rangle =\int\limits_{-1}^1 f(x)g(x)\, dx$. We go up to a cubic because the roots will be the nodes, and we want three nodes. In particular, we want to compute
\begin{align*}
u_1&=v_1\\
v_2&=u_2-\frac{\langle u_2,v_1\rangle}{\langle v_1,v_1\rangle}v_1\\
v_3&=u_3-\frac{\langle u_3,v_1\rangle}{\langle v_1,v_1\rangle}v_1-\frac{\langle u_3,v_2\rangle}{\langle v_2,v_2\rangle}v_2\\
v_4&=u_4-\frac{\langle u_4,v_1\rangle}{\langle v_1,v_1\rangle}v_1-\frac{\langle u_4,v_2\rangle}{\langle v_2,v_2\rangle}v_2-\frac{\langle u_4,v_3\rangle}{\langle v_3,v_3\rangle}v_3
\end{align*} You can find that \begin{align*}
u_1&=1\\
u_2&=x\\
u_3&=x^2-\frac{1}{3}\\
u_4&=x^3-\frac{3}{5}x
\end{align*}
The roots of $u_4$ are $x_{1,2,3}=0,\pm\sqrt{\frac{3}{5}},$ so these are the nodes.
So, our rule is $$\int\limits_{-1}^1f(x)\, dx=c_1f(0)+c_2f\left(\sqrt{\frac{3}{5}}\right)+c_3f\left(-\sqrt{\frac{3}{5}}\right),$$ and all that we have to do is find $c_1,c_2,c_3.$ We impose exactness on $f(x)=1,x,x^2.$ That is, we require that the rule is an equality for these $f$. If we compute what this means, we get that
\begin{align*}
\text{ if }f(x)=1&:\ \ \int\limits_{-1}^1 \, dx=2=c_1+c_2+c_3\\
\text{ if }f(x)=x&:\ \ \int\limits_{-1}^1x \, dx=0=c_2\sqrt{\frac{3}{5}}-c_3\sqrt{\frac{3}{5}}\\
\text{ if }f(x)=x^2&:\ \ \int\limits_{-1}^1x^2 \, dx=\frac{2}{3}=c_2\frac{3}{5}+c_3\frac{3}{5}
\end{align*}
This is a system of linear equations, and you can solve it and find that
\begin{align*}
c_1&=\frac{8}{9}\\
c_2&=\frac{5}{9}\\
c_3&=\frac{5}{9}\\
\end{align*}
This is consistent with the rule that you were looking for.
You can find a proof of the general method of Gaussian quadrature here: https://www.johndcook.com/OrthogonalPolynomials.pdf
Best Answer
Using the scalar product $$\langle f,g\rangle=\sum_{k=1}^5 w(x_k)f(x_k)g(x_k)$$ on the sample space gives the linear system $$ \pmatrix{5&0&2.5\\0&2.5&0\\2.5&0&2.125} \pmatrix{c_0\\c_1\\c_2} = \pmatrix{6\\0\\4.5} $$ leading to $c_1=0$. Eliminating $c_0$ gives $1.75c_2=3\implies c_2=\frac{12}7=1.7142857..$ and lastly $5c_0=6-\frac{30}7=\frac{12}7\implies c_0=\frac{12}{35}=0.342857..$.