Limit and continuity of a multivariable function
The basic concept remains the same. However, with 2D things are more complicated than in 1D. In 1D you either walk forwards or backwards. In 2D we can circle around a point, meaning that, sometimes, the limit may not exist in one direction while it does in another. According to the unit circle we'd have at least 360 different paths to take to reach the point.
For one variable we take one step to the right or one step to the left. For two variables we can take one step up or down, in addition to left and right. For three variables we can take one step to the front or back, in addition to the previous directions. In 2D we have a circle of points around the point which we are evaluating a limit at. In 3D we have a sphere of points.
The idea of 360° for functions of two variables is more or less the "Penrose's Triangle"'s idea. From a certain angle it appears to be a continuous shape. However, from another angle it's revealed that the shape is discontinuous. Use a software to plot a discontinuous function of two variables and spin the viewpoint in an attempt to "hide" the discontinuity.
We have an equation of a circle in 2D and the equation of a sphere in 3D. The equation for the circle is: [math]\displaystyle{ (x - x_0)^2 + (y - y_0)^2 = \delta^2 \iff \delta = \sqrt{(x - x_0)^2 + (y - y_0)^2} }[/math] (we aren't interested in a negative radius, we can disregard the negative root). That equation is also the distance between the circle's origin and [math]\displaystyle{ P_2 }[/math]. Suppose that [math]\displaystyle{ P = (a,b) }[/math] is located anywhere in that circle, excluding the circle's perimeter. Its image, [math]\displaystyle{ f(P) }[/math] is going to be located anywhere in [math]\displaystyle{ \left[L - \epsilon, L + \epsilon\right] }[/math].
Notice how the figure is also a graphical depiction of the property: [math]\displaystyle{ |a - b| = \sqrt{(a - b)^2} }[/math]. Distance cannot be negative. We can view the coordinates of the points as displacement vectors. Both points, [math]\displaystyle{ P_2 }[/math] and the circle's origin, being displaced from the origin of the Cartesian plane to their positions shown in the graph (sum a vector with a point). The radius of that circle can also be interpreted as [math]\displaystyle{ ||\overrightarrow{P_2} - \overrightarrow{C}|| = \delta }[/math], where [math]\displaystyle{ \overrightarrow{C} }[/math] is displacement vector for the circle's origin.
Notation is really the same idea from limits of single variable functions:
[math]\displaystyle{ \lim_{(x,\ y) \ \to \ (x_0,\ y_0)} f(x,y) = L }[/math] (same for any number of variables)
For each [math]\displaystyle{ \epsilon \gt 0 }[/math], there is a [math]\displaystyle{ \delta \gt 0 }[/math], such that every [math]\displaystyle{ (x,\ y) \in D_f }[/math], [math]\displaystyle{ 0 \lt \sqrt{(x - x_0)^2 + (y - y_0)^2} \lt \delta \implies |f(x,y) - L| \lt \epsilon }[/math]. (some textbooks replace the square root with a norm and difference between the coordinates, it's really the same thing)
The concept is virtually the same used for single variable functions. We are considering the smallest distance between two points in 2D that is as close as possible to zero. While the error, the distance between the image and the limit, is the lowest possible value. Note that the definition of a limit for many variables is not considering the path to the point. The concept of one sided limits for many variables is a bit more complicated because when we expand to 2D and 3D there is up and down, front and back, there are way more sides and directions to account for.
Terminology: when a limit does exist at a point [math]\displaystyle{ (x_0, y_0) }[/math], that point is called a limit point or an accumulation point. That term comes from topology.
Continuity: the discussion is exactly the same for one or many variables.
[math]\displaystyle{ \lim_{(x,\ y) \ \to \ (x_0,\ y_0)} f(x,\ y) = f(x_0,\ y_0) }[/math]
If the function is defined at [math]\displaystyle{ (x_0, y_0) }[/math] and the limit converges to that point, the function is continuous at that point. Extend the same reasoning to any arbitrary [math]\displaystyle{ (x, y) }[/math] point chosen from the function's domain and the function is continuous everywhere in its domain.
Squeeze theorem for many variables
It's the same concept from single variable functions, there is no difference.
If [math]\displaystyle{ f(x,\ y) \leq g(x,\ y) \leq h(x,\ y) }[/math] for [math]\displaystyle{ 0 \leq \sqrt{(x - x_0)^2 + (y - y_0)^2} \lt \delta }[/math] and
[math]\displaystyle{ \lim_{(x,\ y) \ \to \ (x_0,\ y_0)} f(x,y) = L = \lim_{(x,\ y) \ \to \ (x_0,\ y_0)} h(x,\ y) }[/math]
Then
[math]\displaystyle{ \lim_{(x,\ y) \ \to \ (x_0,\ y_0)} g(x,\ y) = L }[/math]
The same property of limits for one variable applies to many variables. If a bounded function is multiplied by another with a limit that goes to zero, we can say that the limit of the product is zero too.
"Multi-sided" limits
It's impractical to calculate a limit a hundred times just to check whether it exists or not. We need something to deal with this scenario. It turns out that we have to resort to the knowledge of trajectories, parametric equations. When we plot level curves for functions of two variables we have a path and an equation with two variables, under the special condition that the equation always keep the level a constant. What we are looking for isn't really a vector valued function that describes a trajectory in 3D, the path over the surface of the two variable function. Our problem is that we have a starting point anywhere on the function's domain, the XY plane or part of it, and we wish to walk towards another point, the point that we wish to check whether the limit exists or not. It's hard to picture it, but every step we take in any direction on the function's domain, is reflected on the function's graph.
(Note that we have two straight lines on the domain, parallel to each axis, to reach that point. On the function's surface the trajectories aren't flat though, unless the function itself is a plane)
Which are the easiest paths to check? The XY axis are the most obvious ones. We keep one variable equal to zero and reduce the limit of two variables to a single variable. The next one is to make one variable equal to the other [math]\displaystyle{ x = y }[/math], which translates to walking over the diagonal of the XY plane. Another common and easy path to take is [math]\displaystyle{ y = x^2 }[/math] which is a parabola.
I have a textbook that describes the same concept with a slightly more complicated way. Rather than considering [math]\displaystyle{ x = y }[/math] for example, it considers a curve [math]\displaystyle{ \gamma(t) = (x(t), \ y(t)) }[/math]. It's a function of one variable, in this case time, that outputs positions on the Cartesian plane. That's why we can take [math]\displaystyle{ f(x,y) }[/math] and write [math]\displaystyle{ f(\gamma(t)) }[/math]. We didn't "erase" the two variables, we just "hid" them in another function that traces trajectories on the XY plane. Some examples out there and in some textbooks may have this expression [math]\displaystyle{ f(t,t^2) }[/math]. It's really the same concept of writing [math]\displaystyle{ f(x, y=x^2) }[/math].
From a certain perspective we are using a technique that is looking at the problem of a limit of a function of two variables by imposing certain conditions that allow us to treat it as a limit of a single variable function. When we make [math]\displaystyle{ x = 0 }[/math] or [math]\displaystyle{ y = 0 }[/math] we aren't really transforming a function of two variables into a single variable one. What we are doing is "slicing" the function with a plane and looking at its "silhouette" left by intersecting it with that plane. That silhouette can be interpreted as a single variable function because it's really a line with zero thickness.
I'm going to show some graphs plotted with Geogebra 3D to help get a better view of limits of functions of two variables:
Let's take the classic single variable function [math]\displaystyle{ f(x) = \frac{1}{x} }[/math] and plot it in 3D by keeping the second variable a constant. Note that in this case, the function's behaviour is exactly the same as its single variable counterpart. Approaching from the left or right is going to yield different limits.
Now let's take a look at [math]\displaystyle{ f(x,y) = \frac{x^2}{x^2 + y^2} }[/math]. The blue axis is Z, red is X and green is Y. Notice that if we walk along the X and Y axis the limits at the origin are at a different Z coordinate. If we walk along the diagonals, yet another different Z close to the origin. It's impossible for the function to have different values for the same point in its domain. The graph is rendered as a continuous surface but this is because it's impossible to render at infinite resolution. The discontinuity at the origin has the size of a pixel.
We can very well extend the same reasoning to three and more variables. The textbooks that I know don't have exercises of limits for more than two variables. I suppose it's just an exhausting process to consider even more directions and paths in 3D and beyond.