Find the maximum of $U (x,y) = x^\alpha y^\beta$ subject to $I = px + qy$ 
Let be $U (x,y) = x^\alpha y^\beta$. Find the maximum of the function $U(x,y)$ subject to the equality constraint $I = px + qy$.

I have tried to use the Lagrangian function to find the solution for the problem, with the equation
$$\nabla\mathscr{L}=\vec{0}$$
where $\mathscr{L}$ is the Lagrangian function and $\vec{0}=\pmatrix{0,0}$.
Using this method I have a system of $3$ equations with $3$ variables, but I can't simplify this system:
$$ax^{\alpha-1}y^\beta-p\lambda=0$$
$$\beta y^{\beta-1}x^\alpha-q\lambda=0$$
$$I=px+qx$$

You can certainly simplify that system!
$$\lambda = \frac{\alpha x^{\alpha - 1}y^\beta} p$$
$$\lambda = \frac{\beta y^{\beta - 1}x^\alpha} q$$
Thus
$$\frac{\alpha x^{\alpha - 1}y^\beta} p = \frac{\beta y^{\beta - 1}x^\alpha} q$$
and
$$q\left(\alpha x^{\alpha - 1}y^\beta\right) = p\left(\beta y^{\beta - 1}x^\alpha\right)$$
You can reduce powers:
$$q\alpha \cdot y^\beta = p\beta \cdot xy^{\beta-1}$$
$$q\alpha \cdot y = p\beta \cdot x$$
So now you have two equations:
$$qy\alpha = px\beta$$
$$I = px + qy$$
We then have $$qy = px\frac\beta\alpha$$
$$I = px + px\frac\beta\alpha = x\cdot p\left(1+\frac\beta\alpha\right)$$
So $$x = \frac I {p\left(1+\frac\beta\alpha\right)}$$
Do similar computations for $y$, and you get $$(x, y) = \left(\frac I {p\left(1+\frac\beta\alpha\right)}, \frac I {q\left(1+\frac\alpha\beta\right)}\right)$$
You can then continue with Lagrange multipliers.
