So, I want to be able to derive the Three Point Formula: $$\int_{-h}^h f(x) dx = \omega_0f(x_0) + \omega_1f(x_1) + \omega_2f(x_2)$$ In order to do this, I need to find the coefficients $\omega_{0,1,2}$ and the values $x_{0,1,2}$. I know that you can find the following system of equations: $$ 2h = \omega_{0} + \omega_{1} + \omega_2$$
$$ 0 = \omega_{0} x_{0} + \omega_{1} x_{1}+ \omega_{2} x_2$$
$$ \frac{2}{3} h^3 = \omega_{0} x_{0}^{2} + \omega_{1} x_{1}^{2} + \omega_{2} x_{2}^2 $$
$$ 0 = \omega_{0} x_{0}^{3} + \omega_{1} x_{1}^{3}+ \omega_{2} x_{2}^{3} $$
$$ \frac{2}{5} h^5 = \omega_{0} x_{0}^{4} + \omega_{1} x_{1}^{4}+ \omega_{2} x_{2}^{4} $$
$$0 = \omega_{0} x_{0}^{5} + \omega_{1} x_{1}^{5} + \omega_{2} x_{2}^{5} $$
However, even so far I'm not 100% sure that I have figured out the system of equations correctly. Can anyone provide some guidance for how I'd derive the Gaussian 3-point formula? I have heard that it is useful to use symmetry here but I am unsure how this is applicable. I am also told that $x_1 = 0$, if that helps. Any help would be appreciated.
First, rescale the integral to get $$ I[f] = \frac{1}{h}\int_{-1}^1 f\left(\frac{x}{h}\right)dx\\ \approx \frac{1}{h}\sum_{i=0}^2 w_i f\left(\frac{x_i}{h}\right) $$ in which case the equations to solve for $w_i$ are the first three (scaled) equations you have above $$ \sum_{i=0}^2 w_i = 2\\ \sum_{i=0}^2 w_i x_i = 0\\ \sum_{i=0}^2 w_i x^2_i = \frac{2}{3} $$
This can be recast into a simpler form by using the Legendre polynomials $$ \sum_{i=0}^2 w_i L_j(x_i) = 2\delta_{ij} $$ which is a general form that is true for higher orders of integration, and also for other methods of Gaussian quadrature -- we are doing Gauss-Legendre here.
The weights can then be solved in terms of $\{x_0,x_1,x_2\}$, at which point you are pretty much free to choose the abscissa $x$'s. However, Gaussian quadrature uses the zeros of the $3^\textrm{rd}$ of the associated orthogonal polynomial because it results in a higher accuracy. For $N^\textrm{th}$ order Gaussian quadrature, all polynomials of order $2N-1$ are integrated exactly. The first $N$ by virtue of the fact that you had solved the $N$ equations for the $w_i$ and the higher order by the choice of $x_i$. Since $x_i$ are chosen to be the roots of $L_N(x)$, and any polynomial of degree at most $2N-1$ can be factored into a polynomial of the form $$ p(x) = L_N(x)q(x) + r(x) $$ where the quotient $q(x)$ and remainder $r(x)$ are polynomials of degree at most $N-1$, the first term vanishes and the second term is integrated exactly.
The $3^\textrm{rd}$ order Legendre polynomial is $$ L_3(x) = \frac{1}{2}\left(3x^2-1\right) $$ with roots $x_i\in \{-\sqrt{\frac{3}{5}},0, \sqrt{\frac{3}{5}}\}$. Plugging these into the first three equations leads to the weights $w_i=\{\frac{5}{9},\frac{8}{9},\frac{5}{9}\}$.