This is what I have so far (note this is scratch work):
We have,
$$ \begin{align} \rvert f(x,y,z)-f(1,1,1)\lvert &=\lvert x^2 y+2xz^2-3\lvert\\ &=\lvert x^2y-y+y-1+2xz^2-2x+2x-2\lvert\\ &\leq \lvert y\lvert\lvert x^2-1\lvert+\lvert y-1\lvert+2\lvert zx\lvert \lvert z-1\lvert+2\lvert x-1\lvert\\ &=\lvert y\lvert\lvert x-1\lvert\lvert x+1\lvert+\lvert y-1\lvert+2\lvert zx\lvert \lvert z-1\lvert+2\lvert x-1\lvert\\ &< \lvert y\lvert\lvert x+1\lvert \delta+\:\delta\: + 2\delta\lvert z\lvert\lvert x\lvert+\:2\delta\\ \end{align} $$ This is were I am stuck. I know that since $x-1<\delta$, $y-1<\delta$, and $z-1<\delta$ we can manipulate to have, $x<\delta+1$, $\:y<\delta+1$, $\:z<\delta+1$, $\:x+1<\delta+2.$ But this is just getting messy. I was thinking of defining $\delta$ to be a minimum but would it be a minimum of the deltas above?
Extracting $\delta$, your expression is equal to $$ \delta(|y||x+1| + 2|z||x| + 3) $$ Now, if we say that $\delta$ will at most be chosen to be $1$ (there is nothing special about $1$, any value will do; I just like to use $1$), then we have $$ |y||x+1|\leq 6\\ |z||x| \leq 4 $$ so we get $$ \delta(|y||x+1| + 2|z||x| + 3) \leq \delta(6+8+3) = 17\delta $$ Thus, choosing $\delta = \min\left(\frac{\epsilon}{17}, 1\right)$ will make sure that as long as $$ |x-1|, |y-1|, |z-1| < \delta $$ we have $$ |f(x, y, z) - f(1, 1, 1)| < \epsilon $$
Note, however, that this isn't exactly the definition of continuity. The definition of continuity is that for any $\epsilon>0$ there is a $\delta>0$ such that for any $x, y, z$ with $$ \sqrt{(x-1)^2 + (y-1)^2 + (z-1)^2}<\delta $$ we have $$ |f(x, y, z) - f(1, 1, 1)|<\epsilon $$ You are picking $x, y, z$ from a cube centered around $(1,1,1)$, while the definition wants you to use a ball. In order to apply the above proof to the actual definition you will have to check that this doesn't lead to any kind of problems.