This is a late response to your question however, it never hurts know more math!
In order to calculate $\int_0^1 x^{4}f(x)\,dx$ you use the Method of Undetermined Coefficients. That is:
Let $f_{i}(x)=x^{i}$. Then, calculate $c=\int_0^1 x^{4}f_{i}(x)\,dx$ and $A_{0} f_{i}(x_{0})+A_{1}f_{i}(x_{1})$. From this, you will get an equation of the form $c=A_{0} x_{0}^{i}+A_{1}x_{1}^{i}$. In this problem, we have 4 unknowns hence to solve for $A_{0}, x_{0}, A_{1}, x_{1}$ you must create 4 equations using this method. From there you can solve for your variables.
The Method of Undetermined Coefficients can work with any weight function, using the integral $\int_a^b w(x)f(x)\,dx$ for each $f_{i}(x)=x^{i}$ as defined in the method above.
Let's look first at $\int_{-1}^{1} f(t) \ dt$.
We want the formula to evaluate $\int_{-1}^{1} dt, \int_{-1}^{1} t \ dt, \int_{-1}^{1} t^{2} \ dt, \int_{-1}^{1} t^{3} \ dt, \int_{-1}^{1} t^{4} \ dt, \int_{-1}^{1} t^{5} \ dt, \int_{-1}^{1} t^{6} \ dt $ and $\int_{-1}^{1} t^{7} \ dt$ exactly.
That leads to the following system of equations:
$$ 2 = \omega_{0} + \omega_{1}$$
$$ 0 = \omega_{0} x_{0} + \omega_{1} x_{1}+ \omega_{2}+\omega_{3}$$
$$ \frac{2}{3} = \omega_{0} x_{0}^{2} + \omega_{1} x_{1}^{2} + 2 \omega_{2} x_{2} + 2 \omega_{3} x_{3}$$
$$ 0 = \omega_{0} x_{0}^{3} + \omega_{1} x_{1}^{3}+ 3 \omega_{2} x_{2}^{2} + 3 \omega_{3} x_{3}^{2} $$
$$ \frac{2}{5} = \omega_{0} x_{0}^{4} + \omega_{1} x_{1}^{4}+ 4 \omega_{2} x_{2}^{3} + 4 \omega_{3} x_{3}^{3} $$
$$0 = \omega_{0} x_{0}^{5} + \omega_{1} x_{1}^{5} + 5 \omega_{2} x_{2}^{4} + 5 \omega_{3} x_{3}^{4} $$
$$\frac{2}{7} = \omega_{0} x_{0}^{6} + \omega_{1} x_{1}^{6} + 6 \omega_{2} x_{2}^{5} + 6\omega_{3} x_{3}^{5} $$
$$ 0 = \omega_{0} x_{0}^{7} + \omega_{1} x_{1}^{7} + 7 \omega_{2} x_{2}^{6} + 7\omega_{3} x_{3}^{6}$$
Solve the system using a numerical solver.
Then use the the fact that $$\int_{a}^{b} f(x) \ dx = \frac{b-a}{2} \int_{-1}^{1} f \Big(a+(1+t)\frac{b-a}{2} \Big) \ dt$$
Best Answer
I think you should follow the same procedure you have done above. $$ \begin{align*} \int_{0}^{1}x^{2}f(x)dx & \thickapprox w_{1}f(x_{1})+w_{2}f(x_{2})+w_{3}f(x_{3}) \end{align*} $$ $$ \begin{align*} f(x)=1\Longrightarrow w_{1}+w_{2}+w_{3} & =\int_{0}^{1}1.x^{2}dx=1/3\\ f(x)=x\Longrightarrow w_{1}x_{1}+w_{2}x_{2}+w_{3}x_{3} & =\int_{0}^{1}x.x^{2}dx=1/4\\ & \vdots\\ f(x)=x^{5}\Longrightarrow w_{1}x_{1}^{5}+w_{2}x_{2}^{5}+w_{3}x_{3}^{5} & =\int_{0}^{1}x^{5}.x^{2}dx=1/8 \end{align*} $$ So you have a nonlinear system of equations. You should solve it with a solver. $$ \begin{align*} \int_{a}^{b}w(x)f(x)dx & \thickapprox \sum_{j=1}^{n} w_{j}f(x_{j}) \end{align*}\qquad (*) $$
Additionally as i know in $(*)$ the node points $x_i$ are the zeros of the n-th degree orthogonal polynomial on $[a, b]$ with respect to weight function $w(x)$.
Also after some research on google as i understood there are some packages( for example "ORTHPOL") which can give Gauss quadrature rules for arbitrary weight functions.