What is the most used method for proving continuity for simple functions such as $f(x) = x^{1/3}$
Just for kicks, here is a completely different way to prove that $f(x)=x^2$ is continuous. The identity function is continuous (easy to see, think for example inverse images of open sets). Baby Rudin 4.4 teaches us that if $\lim_{x\to p} f(x)=A$ and $\lim_{x\to p} g(x)=B$, then $\lim_{x\to p} f(x)g(x)=AB$. From this we can infer that if $f$ is continuous at $p$ (meaning $A=f(p)$) and $g$ is continuous at $p$ (meaning $B=g(p)$), then $fg$ is continuous at $p$. Let $f$ and $g$ both be the identity function and conclude that $x^2$ is continuous.
For $x^{1/3}$ we can do the following: $x\to x^3$ is continuous for the same reason that $x\to x^2$ is. Baby Rudin 4.17 tells us that the inverse function of a continuous function is continuous provided the domain is compact. Apply Baby Rudin 4.17 to $x^3$ on the compact sets $[-n,n]$ for all integers $n$, and conclude that $x^{1/3}$, the inverse function of $x^3$, is continuous at all points in $\mathbb{R}$.
Let's look at $f(x) = x^2$ to get the feel for how the proofs work. Fix $\epsilon >0$ and $x \in (-L,L)$, then let's consider $y$ s.t. $| x-y | < \delta $ (you may like to think of delta as a function of $\epsilon$ and $x$ that we have to find). Then we see $$|f(x) - f(y)| = |x^2 - y^2|= | x-y| | x+y| \leq 2L \delta $$ Thus if we take delta like $\delta = \frac{ \epsilon}{2L} $, we see that $$ |f(x) - f(y)| < \epsilon \quad \text{whenever} \quad | x-y | < \delta$$ This applies for all $\epsilon>0$ since we found $\delta(\epsilon)$.$f(x) = x^{1/3}$ is dealt with in the same manner, but instead of having a difference of squares, we have $$ x^3 - y^3 = (x-y ) ( x^2 +xy + y^2 )$$ The trick is to get our delta bound somehow in the difference of $f(x) - f(y)$.
Let's review the $\delta-\epsilon$ proof. We have
$$\begin{align} \left|x^{1/3}-y^{1/3}\right|&=\left|(x^{1/3}-y^{1/3})\left(\frac{x^{2/3}+x^{1/3}y^{1/3}+y^{2/3}}{x^{2/3}+x^{1/3}y^{1/3}+y^{2/3}}\right)\right|\\\\ &=\left|\frac{x-y}{x^{2/3}+x^{1/3}y^{1/3}+y^{2/3}}\right| \tag 1 \end{align}$$
CASE $1$:
If $y=0$, then given $\epsilon>0$, $|x^{1/3}|<\epsilon$ whenever $|x|<\delta=\epsilon^3$.
CASE $2$:
If $y>0$ (we leave the case for which $y<0$ to the reader), then we first choose $\delta<y/2$. Then, $|x-y|<\delta =y/2\implies y/2 < x<3y/2$.
Thus, we have from $(1)$ that given $\epsilon>0$
$$\begin{align} \left|x^{1/3}-y^{1/3}\right|&=\frac{|x-y|}{(y/2)^{2/3}+(y/2)^{1/3}y^{1/3}+y^{2/3}} \\\\ &<\frac{|x-y|}{y^{2/3}}\\\\ &<\epsilon \end{align}$$
whenever $|x-y|<\min\left(y/2,y^{2/3}\epsilon \right)$.
And we are done!