This series illustrate how to derive the linearization from the tangent lines.

We graph f(x, y) = sin(x) - sin(y) and want to find the linearization at (3, 4, sin(3)-sin(4))

Show[Plot3D[Sin[x]-Sin[y],{x, -10, 10}, {y, -10, 10}],ListPointPlot3D[{{3, 4, Sin[3]-Sin[4]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black]]

Here is a zoomed in picture

Show[Plot3D[Sin[x]-Sin[y],{x, 0, 10}, {y, 0, 10}],ListPointPlot3D[{{3, 4, Sin[3]-Sin[4]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black]]

First we graph the tangent lines/linearizations in the x and y directions

Show[Plot3D[Sin[x]-Sin[y],{x, 0, 10}, {y, 0, 10}],ListPointPlot3D[{{3, 4, Sin[3]-Sin[4]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black],ParametricPlot3D[{{3+t, 4, Sin[3]-Sin[4] + t Cos[3]},{3, 4+t, Sin[3]-Sin[4] -t Cos[4]}}, {t, -5, 5}]]

ParametricPlot3D[{{3+t, 4, Sin[3]-Sin[4] + t Cos[3]},{3, 4+t, Sin[3]-Sin[4] -t Cos[4]}}, {t, -5, 5}]

Then we find the plane containing both these lines

Show[Plot3D[{Sin[x]-Sin[y], Sin[3]-Sin[4]+Cos[3](x-3) - Cos[4](y - 4)},{x, -10, 10}, {y, -10, 10}],ListPointPlot3D[{{3, 4, Sin[3]-Sin[4]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black],ParametricPlot3D[{{3+t, 4, Sin[3]-Sin[4] + t Cos[3]},{3, 4+t, Sin[3]-Sin[4] -t Cos[4]}}, {t, -5, 5}]]

Zooming in, we see that it is a better estimate as we get closer

Show[Plot3D[{Sin[x]-Sin[y], Sin[3]-Sin[4]+Cos[3](x-3) - Cos[4](y - 4)},{x, 2, 4}, {y, 3, 5}],ListPointPlot3D[{{3, 4, Sin[3]-Sin[4]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black],ParametricPlot3D[{{3+t, 4, Sin[3]-Sin[4] + t Cos[3]},{3, 4+t, Sin[3]-Sin[4] -t Cos[4]}}, {t, -5, 5}]]

Animate[Show[Plot3D[{Sin[x]-Sin[y], Sin[3]-Sin[4]+Cos[3](x-3) - Cos[4](y - 4)},{x, 3-h,3+h}, {y, 4-h, 4+h}],ListPointPlot3D[{{3, 4, Sin[3]-Sin[4]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black],ParametricPlot3D[{{3+t, 4, Sin[3]-Sin[4] + t Cos[3]},{3, 4+t, Sin[3]-Sin[4] -t Cos[4]}}, {t, -5, 5}]],{h, 10, 0}]

This series illustrates a function where the partials f_x and f_y are well defined at a point, but the linearization remains a poor approximation.

Set f(x,y) =({{{, ({{((x y)/((x^2)+(y^2))), (x,y)≠ (0,0)}, {0, (x, y)=(0,0)}})}})

Plot3D[x y/(x^2+y^2), {x, -5, 5}, {y, -5, 5}]

We can compute the partials
(Subscript[f, x])(0,0)=0
(Subscript[f, y])(0,0)=0
So the linearization is 
L(x, y) = 0

Plot3D[{x y/(x^2+y^2), 0}, {x, -5, 5}, {y, -5, 5}]

Animate[Plot3D[{x y/(x^2+y^2), 0}, {x, -h, h}, {y, -h, h}],{h,10, 0}]

That example was strange becasue f was not continuous at (0,0), so the linearization didn’t make much sense.  
Set f(x, y) = ({{{, ({{((y^3)/((x^2)+(y^2))), (x,y) ≠ (0,0)}, {0, (x, y) = (0,0)}})}})
This function is continuous at (0,0) 
f(r cos θ, r sin θ) = r (sin^3)θ

Plot3D[y^3/(x^2+y^2), {x, -5, 5}, {y, -5, 5}]

And we can compute partials
(Subscript[f, x])(0,0) = 0
(Subscript[f, y])(0,0)=1
So the linearization is
L(x, y) = y

Plot3D[{y^3/(x^2+y^2), y}, {x, -5, 5}, {y, -5, 5}]

Show[Plot3D[{y^3/(x^2+y^2), y}, {x, -5, 5}, {y, -5, 5}],ParametricPlot3D[{{t, 0, 0},{0, t, t}}, {t, -5, 5}]]

Nonetheless, the linearization is never a great approximation.

Animate[Plot3D[{y^3/(x^2+y^2), y}, {x, -h, h}, {y, -h, h}],{h,10, 0}]

I am on a roller coaster with position given by {Sin[t], Cos[3t], Sin[7t]}.  My fear is given by the sum of my height and distance from my security blanket at the origin (z + x^2 + y^2+z^2).  The first cell animates my path along the roller coaster and the second one colors my fear (fear is sometimes negative).

Animate[Show[ParametricPlot3D[{Sin[t], Cos[3t], Sin[7t]}, {t, 0, 2 Pi}], ListPointPlot3D[{{Sin[x], Cos[3x], Sin[7x]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black]], {x, 0, 2Pi}, AnimationRunning->False, AnimationRate->.1]

Animate[Show[ParametricPlot3D[{Sin[t], Cos[3t], Sin[7t]}, {t, 0, 2 Pi}, ColorFunction->Function[{x,y,z},Hue[(z + x^2 + y^2+z^2)/5]]], ListPointPlot3D[{{Sin[x], Cos[3x], Sin[7x]}}, Filling->Bottom, PlotStyle->PointSize[Large], PlotStyle->Black]], {x, 0, 2Pi}, AnimationRate->.1]