Consider a polynomial equation $p(x)=0$ the solutions to which are $x_1,x_2,...,x_n$ (where $p$ is a polynomial function from $\mathbb{C} $ to $\mathbb{C}$). If I wanted to find another polynomial equation q(y)=0 whose solutions are $f(x_1),f(x_2),...,f(x_n)$ for $f:\mathbb{C} \rightarrow \mathbb{C}$, then I would use the substitution $y=f(x) \implies x=f^{-1}(y)$ and substitute this into the original equation which I would then rearrange to get the desired polynomial equation. But what about cases like $y=x^2$ where $f$ is not bijective unless we restrict it's domain? Up until now I have always just written $x=y^{\frac{1}{2}}$ and rearranged to get the desired equation but wouldn't $x=-y^{\frac{1}{2}}$ be just as valid? I tried an exercise of this style using $x=-y^{\frac{1}{2}}$ and got the correct equation.
I would like to understand why both substitutions seem to work the same in getting to the correct final equation. Will this method work generally regardless of the fact that $f$ has different inverses for different parts of this domain. Will this ever cause problems?
Both are valid, and you would need both in general. For a simple example, consider: $$ p(x) = x^2 - x - 2 = (x+1)(x-2) $$
To find the equation(s) with roots $x_1^2, x_2^2$, substitute $y = x^2 \iff x = \pm \sqrt{y}$ with $y \ge 0$. This gives the following two equations to solve.
$$y - \sqrt{y} - 2 = 0 \quad\iff\quad y = 4 \tag{1}$$
$$y + \sqrt{y} - 2 = 0 \quad\iff\quad y = 1 \tag{2}$$
Neither one of the equations has the complete solution set of $\{1,4\}$. To get that, you need to effect both $x = \pm \sqrt{y}$ substitutions, solve the respective equations, then combine the solution sets.
What you did was probably take either one of the equations, move the square root to one side, then square the equation to eliminate the radical, so you ended up with $\,(y-2)^2 = y\,$. This does work, because it is in fact the same equation you'd get by multiplying $(1) \times (2)$:
$$ 0 = \big((y - 2) - \sqrt{y}\big)\big((y - 2) + \sqrt{y}\big) = (y-2)^2 - y $$
You can get the same equation without involving radicals in the process by using polynomial resultants, which amounts to eliminating $x$ between the two algebraic equations $p(x)=0, y=x^2$ in order to obtain an algebraic equation in $y$. For the example here (courtesy WA): $$ \operatorname{res}_x\big(x^2 - x - 2, x^2 - y\big) = y^2 - 5 y + 4 $$
[ EDIT ] More generally, the same approach works with higher-order substitutions $y=x^m$ as well (though not necessarily with arbitrary non-bijective functions). In that case the left-inverse is $x = \omega^k \sqrt[m]{y}$ where $\sqrt[m]{y}$ is the principal value of the $m^{th}$ root, $\omega$ is a primitive $m^{th}$ root of unity, and $k \in \{0, 1, \ldots, m-1\}$. This gives $m$ equations $f\big(\omega^k \sqrt[m]{y}\big) = 0$ for the $m$ different $k$. The union of their solution sets is $\{x_j^m \;\big|\; j = 1, 2, \dots, n\}$, and is the solution set of their product $g(y) = f\big(\sqrt[m]{y}\big) f\big(\omega \sqrt[m]{y}\big) \dots f\big(\omega^{m-1} \sqrt[m]{y}\big)$. From the definition, $g(x)$ is a symmetric polynomial in $\sqrt[m]{y}, \omega \sqrt[m]{y}, \ldots, \omega^{m-1} \sqrt[m]{y}$, so it can be expressed in terms of the elementary symmetric polynomials in $\sqrt[m]{y}, \omega \sqrt[m]{y}, \ldots, \omega^{m-1} \sqrt[m]{y}$, which are the coefficients of $x^m - y$, in fact all $0$ except for the product which is $\pm y$. It follows that $g(y)$ is a real polynomial of degree $n$ which has $\{x_j^m \;\big|\; j = 1, 2, \dots, n\}$ as roots.
For example, using the same $p(x) = x^2 - x - 2$ as above, if $m = 3$ and $\omega$ is a complex cube root of unity then (as verified in WA): $$ \begin{align} g(y) &= \big(\sqrt[3]{y}^2 - \sqrt[3]{y} - 2\big)\big(\omega^2\sqrt[3]{y}^2 - \omega\sqrt[3]{y} - 2\big)\big(\omega\sqrt[3]{y}^2 - \omega^2\sqrt[3]{y} - 2\big) \\ &= y^2 - 7 y - 8 \\ &= (y + 1)(y - 8) \end{align} $$