What are the steps one needs to follow to find a quadrature formula of a certain shape with maximal degree of precision.
For example: Find a quadrature formula of the following shape
$\int_1^2 (x^4 - 1)f(x)dx = A_1f(x_1) + A_2f(x_2) + A_3f(x_3) + R_3(f)$
Should I start by changing the interval to [-1, 1] and then somehow change the formula to a known orthogonal polynomial, Legendre for example ?
You need to construct the Gaussian quadrature for a specific weight function, $w(x) = x^4 - 1$ in your case.
Start with constructing polynomials orthogonal with that weight, i.e. $$ \int_1^2 P_m(x) P_n(x) w(x) = C_n \delta_{nm},\qquad \operatorname{deg} P_n(x) = n. $$ Since you need a quadrature with three points, the nodes of the quadrature are roots of the $P_3(x)$ polynomial. The easiest way to build the orthogonal set by hand is to apply Gramm-Schmidt process to the $1, x, x^2, x^3$ set.
$P_0(x)$ is constant and could be taken as $1$ without loosing any generality. $P_1(x) = x - \alpha$. Since $P_1(x)$ should be orthogonal to $P_0(x)$ we have $$ \int_1^2 P_1(x) P_0(x) (x^4-1)dx = \int_1^2 (x-\alpha)(x^4-1) dx = 9 - \frac{26}{5}\alpha = 0, $$ so $\alpha = \frac{45}{26}$.
Next, $P_2(x) = x^2 - \beta x - \gamma$. Integrating with $P_0(x)$ and $P_1(x)$ we have $$ \int_1^2 P_2(x) P_0(x) (x^4-1)dx = \frac{332}{21} - 9\beta -\frac{26}{5}\gamma = 0\\ \int_1^2 P_2(x) P_1(x) (x^4-1)dx = \frac{555}{728} - \frac{127\beta}{546} = 0, $$ so $\beta = \frac{1665}{508}, \gamma = -\frac{56165}{21336}$. (Ewww)
Finally, $P_3(x) = x^3 - \delta x^2 - \epsilon x - \zeta$. $$ \int_1^2 P_3(x) P_0(x)(x^4-1)dx = -\frac{332 \delta }{21}-9 \varepsilon -\frac{26 \zeta }{5}+\frac{225}{8}\\ \int_1^2 P_3(x) P_1(x)(x^4-1)dx = -\frac{555 \delta }{728}-\frac{127 \varepsilon }{546}+\frac{17783}{9360}\\ \int_1^2 P_3(x) P_2(x)(x^4-1)dx = \frac{18419}{284480}-\frac{40337 \delta }{2987040} $$ And the $P_3(x)$ is $$P_3(x) = x^3 - \frac{386799}{80674}x^2 +\frac{4566058}{605055}x - \frac{2512489}{645392}$$ with roots $$ x_1 \approx 1.24311506656362394748261056490...\\ x_2 \approx 1.62867605333517443866228087173...\\ x_3 \approx 1.92280193362526141007198140964... $$
Now we are able to compute the weights by integrating Lagrange basis polynomials $$ w_i = \int_1^2 \ell_i(x) (x^4 - 1) dx, \qquad \ell_i = \prod_{j \neq i} \frac{x-x_j}{x_i - x_j}\\ w_1 \approx 0.498562063242571816059749610...\\ w_2 \approx 2.242930705042292632885672141...\\ w_3 \approx 2.458507231715135551054578248... $$
Checking that it actually works for integrating every polynomial of degree up to five including: $$ \int_1^2 1 (x^4 - 1)dx = 5.2 = w_1 + w_2 + w_3\\ \int_1^2 x (x^4 - 1)dx = 9 = w_1 x_1 + w_2 x_2 + w_3 x_3\\ \int_1^2 x^2 (x^4 - 1)dx = 15.809523 = w_1 x_1^2 + w_2 x_2^2 + w_3 x_3^2\\ \int_1^2 x^3 (x^4 - 1)dx = 28.125 = w_1 x_1^3 + w_2 x_2^3 + w_3 x_3^3\\ \int_1^2 x^4 (x^4 - 1)dx = 50.577777 = w_1 x_1^4 + w_2 x_2^4 + w_3 x_3^4\\ \int_1^2 x^5 (x^4 - 1)dx = 91.8 = w_1 x_1^5 + w_2 x_2^5 + w_3 x_3^5\\ \int_1^2 x^6 (x^4 - 1)dx = 167.9480519 \neq 167.947233 = w_1 x_1^6 + w_2 x_2^6 + w_3 x_3^6. $$