We use the $\delta - \epsilon$ definition, where $\forall \epsilon$, $\exists \delta$ such that $|x-x_0|\lt\delta $ implies $|f(x)-f(x_0)| \lt \epsilon$.
Here is my attempt:
$|x^3-x_0^3| \lt \epsilon$ implies $|x-x_0||x^2+x_0x+x_0^2| \lt \epsilon$. hence, $|x-x_0| \lt \frac{\epsilon}{|x^2-x_0x+x_0^2|}$.
Let $\delta=1$. $|x-x_0| \lt 1$. Hence, $|x| \lt 1+|x_0|$.
From this, I replaced all $x$ in $\frac{\epsilon}{|x^2-x_0x+x_0^2|}$ with $1+|x_0|$. So, $\delta=min(1,\frac{\epsilon}{3x_0^2+3x_0+1})$
Would this be fine? It seems someone posted a question like this but he got $\delta=min(1,\frac{\epsilon}{2x_0^2+3x_0+1})$ with a different method...
It may not seem like it but I'm following the same idea of your proof that it is using the maximum and minimum values of the function. But you have to get the maximum and the minimum of the restricted function at closed intervals. I'm just rewriting your proof in a leaner way.
The case $x_0=0$ is trivial. Without loss of generality, you can supose $x_0>0$ and restrict the $ x $ variable to the interval $$ I=[\frac{1}{2}x_0,\frac{3}{2}x_0]=[x_0-\frac{1}{2}x_0,x_0+\frac{1}{2}x_0 ] $$ The function $x\mapsto (x^2-x_0x+x_0^2)$ is growing in the interval $I$. It have a global minimum in $x=\frac{1}{2}x_0$ and a global maximum in $x=\frac{3}{2}x_0$. Then $$ 0<\frac{3}{4}(x_0)^2= \min_{x\in I} \;(x^2-x_0x+x_0^2)\leq (x^2-x_0x+x_0^2)\leq \max_{x\in I} \;(x^2-x_0x+x_0^2)=\frac{7}{4}(x_0)^2 $$ This means that the expression $ x^2-x_0x + x_0^2 >0$ for all $ x \in I $. Let $M=\max_{x\in I} \;(x^2-x_0x+x_0^2)$. For all $\epsilon>0$ there is $\delta =\epsilon/M>0$ such that \begin{align} 0<|x-x_0|<\delta \implies & 0<|x-x_0|<\epsilon/M \\ \implies & 0<|x-x_0|\cdot M<\epsilon \\ \implies & 0<|x-x_0|\cdot (x^2-x_0x+x_0^2)<\epsilon \\ \implies & 0<|(x-x_0)\cdot (x^2-x_0x+x_0^2)|<\epsilon \\ \implies & 0<|x^3-x_0^3|<\epsilon \\ \end{align}