How does one establish continuity at a point for a multivariable scalar field? I am trying to understand continuity for multivariable scalar fields, $f: \Bbb R^n \rightarrow \Bbb R$. 
In single variable calculus, I understand the $\delta-\epsilon$ formulation 
of continuity.

For a function $f: \Bbb R \rightarrow \Bbb R$, $f$ is continuous if
  for any arbitrary $\epsilon >0$, $\exists \delta$ such that
$$|x - x_0 | < \delta \Rightarrow |f(x) - f(x_0)| <\epsilon$$

For example, suppose $f(x) = 3x - 6$. 
\begin{align*}
|f(x)-f(a)|&<\epsilon\\
|(3x-6)-(3a-6)|&<\epsilon\\
|3x-3a|&<\epsilon\\
3|x-a|&<\epsilon\\
|x-a|&<\frac{\epsilon}{3}\\
\end{align*}
So for $a \in \Bbb R$, then $\delta = \frac{\epsilon}{3}$ implies $|f(x) - f(a)|<\epsilon$
For scalar fields, Tom Apostol says 

A function $f: \Bbb R^n \rightarrow \Bbb R$ is continuous at $\mathbf{a}$ if 
  
  
*
  
*$f$ is defined at $\mathbf{a}$ 
  
*$\displaystyle \lim_{\mathbf{x}\rightarrow \mathbf{a}} f(\mathbf{x})=f(\mathbf{a})$
  

But I don't understand how to apply the $\delta-\epsilon$ formulation here. 
Suppose I want to show $f(x,y) = x^4 + y^4 - 4 x^2 y^2$ is continuous at some point, say $\mathbf{a} = (3,4)$. 
\begin{align*}
|f(x,y)-f(3,4)|&<\epsilon \\
|x^4 + y^4 - 4 x^2 y^2 - 3^4 - 4^4 + 4\times 3^2 4^2|&<\epsilon\\
\end{align*}
But this expression is so much uglier than the single variable case. Should I be approaching this just like the single variable case? In other words, should the final simplification look like
\begin{align*}
||\mathbf{x}-\mathbf{a}|| &<\delta\\
\sqrt{(x-x_0)^2 + (y-y_0)^2}&<\delta\\
\end{align*}
Or is this the wrong approach? I have googled a lot, but I can't find examples that show continuity for multivariable functions, just ones that disprove it. 

Not that bad. For
\begin{align*}
&|f(x,y)-f(3,4)|\\
&=|x^{4}+y^{4}-4x^{2}y^{2}-3^{4}-4^{4}+4\cdot 3^{2}4^{2}|\\
&=|(x^{2}+3^{2})(x+3)(x-3)+(y^{2}+4^{2})(y+4)(y-4)-4(xy+3\cdot 4)(xy-3\cdot 4)|\\
&\leq|x^{2}+3^{2}||x+3||x-3|+|y^{2}+4^{2}||y+4||y-4|+4|xy+3\cdot 4||(x-3)y+3(y-4)|\\
&\leq(x^{2}+9)(|x|+3)|x-3|+(y^{2}+16)(|y|+4)|y-4|+4(|xy|+12)(|x-3||y|+3|y-4|).
\end{align*}
For $(x-3)^{2}+(y-4)^{2}<\delta^{2}=\min\{1,\epsilon\}$, then $|x-3|<1$ and $|y-4|<1$ and we have $|x|+3\leq|x-3|+6<7$, $x^{2}+9<25$, $|y|+4\leq|y-4|+8<9$, $y^{2}+16<41$, $|xy|+12<32$, so 
\begin{align*}
|f(x,y)-f(3,4)|\leq M(|x-3|+|y-4|)\leq 2M\sqrt{|x-3|^{2}+|y-4|^{2}}<2M\epsilon,
\end{align*}
where $M>0$ is a universal constant, essentially, it is the sum of all those bounds, one may take $M=25\cdot 7+41\cdot 9+4\cdot 32(5+3)$.
