Derive the parameters of Guass quadrature with three points $$ \int_{-1}^1 f(x)\,\mathrm dx \approx C_1 f(\xi_1)+C_2 f(\xi_2) + C_3 f(\xi_3) $$ such that the integral is exact up to $x^5$.
Attempt
Since we are looking at three points, then we may choose $\xi_1 = -1 + \frac{2}{4} = - \frac{1}{2}$, $\xi_2 = - \frac{1}{2} + 1/2 = 0 $ and $\xi_3 = \frac{1}{2} $. So we need to find the weights $C_i$. Since we want exact up to $x^5$, we want first
$$ C_1 + C_2 + C_3 = \int_{-1}^1 1\,\mathrm dx = 2 $$
$$ C_1 (-1/2) + C_2 \cdot 0 + C_3 \cdot 1/2 = \int_{-1}^1 x \,\mathrm dx =1 $$
but then if we keep going until the 5th order, we will have 5 equations and three unknowns. Is this the way to solve the problem or am I doing something wrong?
With 3 points and 3 weights, you have 6 unknowns allowing you to satisfy 6 equations. There are 6 monomials $x^j$ for $0\leq j\leq 5$, so we should be able to solve for the weights and nodes that yield the exact values for the integrals for polynomials up to the fifth degree.
The standard way to do this is by using orthogonal polynomials. This is explained in detail here. The idea is simple, if you demand that the quadrature rule is exactly valid for the Legendre polynomials $P_j(x)$ for $0\leq j\leq 3$ and for $x P_3(x)$ and $x^2 P_3(x)$, then because of orthogonality the integrals except for $j = 0$ are zero. That $f(x)$ to be $P_3(x)$, $x P_3(x)$ and $x^2 P_3(x)$ must be zero, implies that unless the nodes must be the zeros of $P_3(x)$, and the weights then follow from the other 3 equations. You don't actually need to solve these equations as there is an equation in terms of $P_3(x)$ and $P_2(x)$ for the weights.
There exists another far more powerful method to compute quadrature rules, which will also work in more complex cases where the formalism using orthogonal polynomials doesn't work. This involves solving the equations you get when you write down that the quadrature rule must work for the first 6 monomials using a suitable generating function.
In this case, we can exploit that odd functions integrate to zero and consider the integrals:
$$\int_{-1}^{1}x^{2j}dx = \frac{2}{2j+1}$$
The quadrature equations are then:
$$\sum_{k=1}^{2}C_k \xi_k^{2j} = \frac{2}{2j+1}$$
for $0\leq j\leq 2$. We impose the extra condition that the nodes are positive as we're only considering even functions, so there will only be 2 nodes one of which will be zero. At the end of the computations, we can add negative nodes and split up the weights symmetrically, this then doesn't affect even functions while odd functions will yield zero.
Multiplying the quadrature equation by $1/z^{j+1}$ and summing over $j$ from zero to infinity (disregarding that the equation is not valid for $j>2$) yields for the left hand side:
$$\sum_{k=1}^{2}C_k \frac{1}{z} \sum_{j=0}^{\infty}\frac{\xi_k^{2}}{z^j} = \sum_{k=1}^{2}\frac{C_k}{z-\xi_k^2} \tag{1}$$
The right hand side is a series expansion in powers of $1/z$, so if we can rewrite that series as a rational function we could get to the above form using partial fractions expansions and extract the nodes and the weights. The Padé approximant method is an efficient way to do exactly that. In this case we have the series:
$$\frac{2}{z} + \frac{2}{3 z^2} + \frac{2}{5 z^3} + \cdots$$
If we can write down a rational function whose series expansion agrees with the first 3 terms of this series that has a partial fraction expansion of the form (1), then we're done. It's convenient to put $u = 1/z$, and work with the series:
$$2 u + \frac{2 u^2}{3} + \frac{2 u^3}{5} + \cdots$$
To simplify things, we divide this by $u$ and find the Padé approximant of that first. We thus seek a rational function that yields the first 3 terms of the series:
$$2 + \frac{2 u}{3} + \frac{2 u^2}{5} + \cdots \tag{2}$$
Dividing by $u$ will have eliminated the node at zero, this means that the degree of the denominator will be $1$. We thus seek a the rational function of the form:
$$\frac{p + q u}{1 + a u}$$
Multiplying the series (2) by the denominator, working out the result to order $u^2$ and equating the result to the numerator, yields:
$$p + q u = 2 +\left(2a + \frac{2}{3}\right) u + \left(\frac{2 a}{3}+\frac{2 u^2}{5}\right) u^2 $$
Equating the coefficient of $u^2$ on the right hand side to zero yields $a = -\dfrac{3}{5}$. We see that $p = 2$, and $q =-\dfrac{8}{15}$. Multiplying the rational function by $u$ and substituting $z = 1/u$ then yields the Padé approximant for the original series as:
$$\frac{2}{z}\frac{z - \frac{4}{15}}{z - \frac{3}{5}}$$
The partial fraction expansion is given by:
$$\frac{8}{9}\frac{1}{z} + \frac{10}{9}\frac{1}{z - \frac{3}{5}}$$
So, we see that the node at zero has a weight of $\dfrac{8}{9}$ while the node at $\sqrt{\dfrac{3}{5}}$ has a weight of $\dfrac{10}{9}$. To revert to the original problem where the function is not constrained to be even, we split up the nodes symmetrically about zero. This then yields nodes at $\pm\sqrt{\dfrac{3}{5}}$ with weights $\dfrac{5}{9}$ and a node at zero with weight $\dfrac{8}{9}$.