Code to Illuminate Multivariable Calculus Topics

Al Hibbard

Central College

Overview

NumericPlot3D

CenteredNumericPlot3D

CenteredContourPlot

Graphics3DViewer

CenteredPlot3D

Plot3DExt

RainbowParametricPlot

LineIntegralIllustration

ShowApproximations3D

NumericPlot3D

g[x_, y_] :=  (x^2) + 3  (y^2)

ContourPlot[g[x, y], {x, -1, 1}, {y, -1, 1}];

NumericPlot3D[g[x, y], {x, -1, 1}, {y, -1, 1}]

NumericPlot3D[g[x, y], {x, -1, 1, 0.5}, {y, -1, 1}]

NumericPlot3D[
      g[x, y], {x, -1, 1, 0.5}, {y, -1, 1, 0.25}]

data = 
      NumericPlot3D[g[x, y], {x, -1, 1, 0.5}, {y, -1, 1, 0.5},  
        Frame -> False]

First[data][[2]] // Rest

Code

CenteredNumericPlot3D

g[x_, y_] :=  (x^2) + 3  (y^2)

CenteredNumericPlot3D[g[x, y], {x, 0}, {y, 0}, 
      Radii -> 1]

NumericPlot3D[g[x, y], {x, -10, 10, 0.5}, {y, -1, 1, 0.5}, 
      Focus -> {0, 0}, Radii -> 1.0]

CenteredNumericPlot3D[g[x, y], {x, 0}, {y, 0}, 
      Radii -> {1, 4.0}]

NumericPlot3D[g[x, y], {x, -10, 10, 0.5}, {y, -1, 1, 0.5}, 
      Focus -> {0, 0}, Radii -> {1.0, 4}]

Table[
      CenteredNumericPlot3D[g[x, y], {x, 1}, {y, -2}, 
        Radii -> (0.5^k)], {k, 1, 5}]

CenteredNumericPlot3D[(x^2) + 3 (y^2), {x, 0},  {y, 0}, 
      Radii -> {(1/2), 1}, Zoomable -> True]

code

CenteredContourPlot

f[x_, y_] := 2 x y  -  Cos[x]  +  Sin[ x  -  y]

ContourPlot[f[x, y], {x, -1, 3}, {y, 0, 4}];

CenteredContourPlot[f[x, y], {x, 1}, {y, 2}, 
        Radii -> 2];

Options[CenteredContourPlot]

CenteredContourPlot[f[x, y], {x, 1}, {y, 2}, 
        Radii -> {2,  1}];

Table[
        ContourPlot[f[x, y], {x, -2, 3}, {y, -4, 5}, 
          Focus -> {2, 3}, Radii -> (2^(-k))], {k, -2, 
          3}];

CenteredContourPlot[f[x, y], {x, 1}, {y, 2}, Radii -> 2, 
      Zoomable -> True];

CenteredContourPlot[2 x y - Cos[x] + Sin[x - y], {x, 1}, {y, 2}, 
        Radii -> 3,  ShowGrid -> False, 
        ColorFunction -> Automatic];

CenteredContourPlot[2 x y - Cos[x] + Sin[x - y], {x, 1}, {y, 2}, 
      Radii -> 3, Zoomable -> True,  ShowGrid -> False, 
      ColorFunction -> Automatic];

ContourPlot[f[x, y], {x, -2, 3}, {y, -4, 5}, 
      Focus -> {2, 3}, Radii -> {3, 1}, 
      Zoomable -> True];

code

Graphics3DViewer

f[x_, y_] := Sin[x + Cos[y]]

gr = 
        Plot3D[f[x, y], {x, -4, 5}, {y, -6, 7}, 
          ColorFunction -> Hue, PlotPoints -> 35];

Graphics3DViewer[gr]

g[x_, y_] := 5 - 3  (x^2) - 2  (y^2)

gr = 
        Plot3D[g[x, y], {x, -6, 6}, {y, -3, 3}, 
          PlotRange -> All];

Graphics3DViewer[gr, Viewer -> "Browser"]

code

CenteredPlot3D

Clear[f]
f[x_, y_] := x y Exp[-((x^2) + (y^2))]

CenteredPlot3D[f[x, y], {x, 1.2}, {y, 1.4}, Radii -> {2, 3},  
        PlotRange -> All];

CenteredPlot3D[f[x, y], {x, 1.2}, {y, 1.4}, Radii -> 2,  
        PlotRange -> All];

Table[
        CenteredPlot3D[f[x, y], {x, 1.2}, {y, 1.4}, 
          Radii -> (2^(-k))], {k, 0, 5}];

CenteredPlot3D[(E^(-(x^2) - (y^2))) x y, {x, 
        1.2},  {y, 1.4}, Radii -> 1, Zoomable -> True, 
      AxesLabel -> {"x", "y", "z"}, 
      PlotRange -> All];

code

Plot3DExt

Clear[f]
f[x_, y_] := x y Exp[-((x^2) + (y^2))]

gr = 
        Plot3DExt[f[x, y], {x, -3, 3}, {y, -3, 3}, 
          Focus -> {1.2, 1.4}, PlotRange -> All];

Graphics3DViewer[gr]

g[x_, y_] := 3  (x^2) - 2  (y^2)

plane = TangentPlane[g[x, y], {x, 2}, {y, -2}]

Plot3DExt[g[x, y], {x, -2, 3}, {y, -3, 3}, 
        PlaneEquation -> plane];

gr = 
        Plot3DExt[g[x, y], {x, -2, 3}, {y, -3, 3}, 
          PlaneColor -> Red, PlaneEquation -> plane];

Graphics3DViewer[gr]

code

RainbowParametricPlot

ParametricPlot[{Cos[3  t], Sin[5  t]}, {t, 0, 2  π},  
        PlotStyle -> RGBColor[0, 0, 1],  
        AspectRatio -> Automatic];

TraceParametricPlot[{Cos[3  t], Sin[5  t]}, {t, 0, 2  π, 8},  
        PlotStyle -> RGBColor[0, 0, 1],  
        AspectRatio -> Automatic];

RainbowParametricPlot[{Cos[3  t], Sin[5  t]}, {t, 0, 2  π}, 16,  
        PlotStyle -> RGBColor[0, 0, 1],  
        AspectRatio -> Automatic];

RainbowParametricPlot[{Cos[3  t], Sin[5  t]}, {t, 0, 2  π}, 16,  
        PlotStyle -> RGBColor[0, 0, 1],  
        AspectRatio -> Automatic,  
        IncludeDots -> True];

code

LineIntegralIllustration

Clear[P, Q]
P[x_, y_] := x - y
Q[x_, y_] := x + y

gr1 = 
        PlotVectorField[{P[x, y], Q[x, y]}, {x, -3, 3}, {y, -3, 3}, 
          ColorFunction -> Hue,  
          DisplayFunction -> Identity];
gr2 = 
        Plot[(Sqrt[1 - (x^2)]), {x, -1, 1}, 
          PlotStyle -> RGBColor[0, 0, 0],  
          DisplayFunction -> Identity];
Show[{gr2, gr1}, AspectRatio -> Automatic,  
        DisplayFunction -> $DisplayFunction];

LineIntegralIllustration[{P[x, y], Q[x, y]}, {Cos[t], Sin[t]}, {t, 
          0,  Pi}, 4,  AspectRatio -> Automatic];

LineIntegralIllustration[{P[x, y], Q[x, y]}, {Cos[t], Sin[t]}, {t, 
          0,  Pi}, 12,  AspectRatio -> Automatic];

code

ShowApproximations3D

Clear[f, g]
f[x_, y_] := 
      x y Exp[-((x^2) + (y^2))]
g[x_, y_] := Sin[x + Cos[y]]

gr = 
        ShowApproximations3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
          PlotRange -> All];

Graphics3DViewer[gr]

Options[ShowApproximations3D]

gr = 
      ShowApproximations3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
        Method -> "closestCorner", Divisions -> {12, 6}, 
        PlotRange -> All, DisplayVolume -> True]

gr = 
        ShowApproximations3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
          Method -> "closestCorner", Divisions -> {12, 6}, 
          PlotRange -> All];

Graphics3DViewer[gr]

Plot3D[g[x, y], {x, -3, 3}, {y, -3, 3}];

gr = 
        ShowApproximations3D[g[x, y], {x, -3, 3}, {y, -3, 3}, 
          Method -> "lowerLeft", Divisions -> {10, 20}, 
          Color -> Yellow];

Graphics3DViewer[gr]

RiemannSum[g[x, y], {x, -3, 3}, {y, -3, 3}, 
      Divisions -> {10, 20}, Method -> "lowerLeft"]

NRiemannSum[g[x, y], {x, -3, 3}, {y, -3, 3}, 
      Divisions -> {10, 20}, Method -> "lowerLeft"]

NRiemannSum[g[x, y], {x, -3, 3}, {y, -3, 3}, 
      Divisions -> {10, 20}, Method -> "midpoint"]

NIntegrate[g[x, y], {x, -3, 3}, {y, -3, 3}]

NRiemannSum[g[x, y], {x, -3, 3}, {y, -3, 3}, 
      Divisions -> {200, 200}, Method -> "lowerLeft"]

g[x_, y_] := (x^2) + (y^2)

\(TableForm[
      Table[{3\^k, NIntegrate[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "lowerLeft"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "lowerRight"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "upperLeft"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "upperRight"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "midpoint"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, 
            Method -> "closestCorner"]}, {k, 0, 4}], 
      TableHeadings -> {None, {"n", "actual", "lowLft", "lowRgt", "upLft", "upRgt", "mid", "clsCrn
"}}]\)

g[x_, y_] := Cos[x  +  y]

\(TableForm[
      Table[{3\^k, NIntegrate[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "lowerLeft"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "lowerRight"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "upperLeft"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "upperRight"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, Method -> "midpoint"], 
          NRiemannSum[g[x, y], {x, \(-3\), 3}, {y, \(-3\), 3}, 
            Divisions -> {3\^k, 3\^k}, 
            Method -> "closestCorner"]}, {k, 0, 4}], 
      TableHeadings -> {None, {"k", "actual", "lowLft", "lowRgt", "upLft", "upRgt", "mid", "clsCrn
"}}]\)

Code to initialize at outset

Preliminaries

Off[General::"spell1"];
Off[
        General::"spell"];
Needs["Graphics`Legend`"]
  
Needs["Graphics`Arrow`"]
Needs["Graphics`Colors`"]
Needs["Graphics`PlotField`"]
Needs["Graphics`FilledPlot`"]
Needs["Utilities`FilterOptions`"]
Needs["JLink`"]

\(SetOptions[Plot, \ 
          PlotStyle\  -> \ {{Red}, {Blue}, {Magenta}, {Cyan}, 
	{Green}, {Yellow}, \ {Black}}];\)
SetOptions[ListPlot,  PlotStyle  ->  {Blue}];
vp  =  {1.429,  .645, 3};
SetOptions[
        Plot3D,  ViewPoint -> vp ];
SetOptions[Graphics3D,  
        ViewPoint -> vp ];
SetOptions[ParametricPlot3D,  
        ViewPoint -> vp ];
Unprotect[
        Black];
Black  =  
        RGBColor[0, 0, 0];
Protect[
        Black];

Code sections

LiveJar

$CurrentNotebookName := "WindowTitle" /. 
        NotebookInformation[EvaluationNotebook[]]
$CurrentNotebookPath := ("FileName" /. 
              NotebookInformation[EvaluationNotebook[]])[[1]] // 
        ToFileName

$JarFileName  =  "live.jar";

    
FileInDirectoryQ[fn_, dir_] := 
      Module[{path,  files,  cwd}, 
cwd  =  
          Directory[]; 
path  =  
          ToFileName[dir]; 
SetDirectory[
          path]; 
files  =  
          FileNames[]; 
SetDirectory[
          cwd]; 
MemberQ[files, 
          fn]]

    
LiveInCurrentDirectoryQ[] := 
      FileInDirectoryQ[$JarFileName, $CurrentNotebookPath]

    
LiveInSomeMMAPathQ[] :=  
      Apply[Or, 
        Map[FileInDirectoryQ[$JarFileName, #] &, $Path]]

    
$LivePresent  =  
        LiveInCurrentDirectoryQ[]  ||  
          LiveInSomeMMAPathQ[];

    
LiveMMAPath[] :=  Select[
            Cases[Map[{#, FileInDirectoryQ[$JarFileName, #]} &, $Path], {_, 
                True}], StringLength[First[#]] > 3 &] // First // 
        First

    
LiveJarPath  :=  
      If[LiveInCurrentDirectoryQ[], $CurrentNotebookPath, 
        If[LiveInSomeMMAPathQ[],  LiveMMAPath[], "none"]]

\(QueryForLivePath[] := 
      Module[{nb}, 
nb\  = \ 
          NotebookCreate[]; 
NotebookWrite[nb, \ 
          Cell["You apparently do not have the 'live.jar' file installed in the same place as the notebook in which you are working nor in any of the standard locations within the Mathematica layouts. As one last chance (before you go to the web and download it - see below, if needed), you can search for it among your directories. When ready, click on the 'Find it' button. Note: If you don't know where it is or you can't find it, just pick any file.", \ "Text"]]; 
$FBPath\  = \ ""; 
NotebookWrite[nb, \ 
          Cell[BoxData[
              ButtonBox["Find live.jar", 
                ButtonFunction :> \((FindLive[#] &)\), 
                ButtonEvaluator -> Automatic]], "Input", 
            Active -> True]];
]\)

\(FindLive[x_] := 
      Module[{ljpath, \ pos}, 
LinkWrite[$ParentLink, \ 
          FrontEnd`FileBrowse[False]]; 
$FBPath\  = 
          LinkRead[$ParentLink]; 
pos\  = 
          StringPosition[$FBPath, "live.jar"]; 
If[
          Length[pos]\  == 0, 
NotebookWrite[
            ButtonNotebook[], \ 
            Cell[TextData[{"You did not found the file. To see live 3D graphics, go to

", \ 
                  ButtonBox["http://www.vis.uni-stuttgart.de/~kraus/LiveGraphics3D/download.html", 
                    ButtonData :> {URL["http://www.vis.uni-stuttgart.de/~kraus/LiveGraphics3D/download.html"], None}, 
                    ButtonStyle -> "Hyperlink"], \ "

to download it. Then place it in the same folder as the file in which you were trying to work. You may now close this window."}], \ "Text"]]; 
          LiveJarPath = "none", 
$FBPath\  = \ 
            StringDrop[$FBPath, 
              First[pos]]; 
LiveJarPath = $FBPath; 
NotebookWrite[ButtonNotebook[], \ 
            Cell[TextData[{"You found the file at 

", $FBPath, "

Great. So you can now view any 3D output of any function by 'feeding' it into the Graphics3DView function. Close this window and begin."}], \ "Text"]]];
]\)

\(OneTimeJLinkStuff := 
      Module[{path}, 
$LiveJarReady\  = \ 
          False; 
path = 
          LiveJarPath; 
LJ = 0; 
If[
          path\  == "none", QueryForLivePath[]]; 
If[
          path\  =!= "none", $LiveJarReady = 
            True; 
InstallJave[]; 
AddToClassPath[LiveJarPath <> $JarFileName]]]\)

\(MMAViewer[graphics_, animOpts___] := 
      AppletViewer["Live", {"INPUT=" <> 
            StringReplace[
              ToString[
                LiveForm[graphics, 
                  animOpts]], """ -> "''"], "WIDTH=600", "HEIGHT=600"}]\)

\(\(\(\ \)\(\(Options[
          Graphics3DViewer] = {Viewer -> "Mathematica"};\)\(\ \)\(
\)
    \( (*\ \(Options[
            Graphics3DViewer] = {Viewer -> "Browser"};\)\ *) \)\)\)\)

\(Graphics3DViewer[graphics_, animOpts___] := 
      Module[{ver, \ 
          viewer}, 
ver\  = \ $VersionNumber; 
viewer = 
          Viewer /. Flatten[{animOpts, 
                Options[Graphics3DViewer]}]; 
 (*\ \(If[
              animOpts =!= Null, 
animOpts = 
                DeleteCases[animOpts, \ 
                  Viewer\  -> \ _]];\)\ *) 
If[$LiveJarReady =!= True, 
          OneTimeJLinkStuff]; 
If[ver ≥ \ 5, 
          If[viewer === "Browser", MMAToBrowser[graphics, animOpts], 
            MMAViewer[graphics, animOpts]], 
          MMAToBrowser[graphics, animOpts]]]\)

\(\(LiveForm[graphics_, animOpts___] := 
        Module[{graphics3d = 
              Switch[graphics, 
_ContourGraphics, 
                  Graphics3D[SurfaceGraphics[graphics]], 
_DensityGraphics, 
                  Graphics3D[SurfaceGraphics[graphics]], 
_SurfaceGraphics, 
                  Graphics3D[graphics], 
_, 
                  graphics] //. SequenceForm[
                    x___] :> 
                  StringForm[StringJoin[Table["``", {Length[{x}]}]], x]}, 
          Switch[graphics3d, 
_List, \(HoldForm[ShowAnimation]\)[
              LiveForm /@ graphics3d, 
              InputForm[N[Flatten[{animOpts}]]]], 
_, 
            InputForm[N[graphics3d]]]];\)\)

MMAToBrowser[graphics_, animOpts___] := 
        Module[{ps, x,  nfilename , mfilename,  
            filename = "LiveJarGraphic" <> 
                ToString[LJ++]}, 
mfilename = 
            filename <> ".m"; 
If[Not[OSChosenQ], 
            QueryOS[]]; nfilename  =  StringJoin[LiveJarPath, mfilename]; 
          ps = Unprotect[Real];  
          Format[x_Real, InputForm] := 
            OutputForm[
              NumberForm[x, 5, 
                NumberFormat -> (If[#3 == "", #1, 
                        SequenceForm[#1, "*^", #3]] &)]];  
          Protect @@ ps; 
WriteString[nfilename, 
            ToString[LiveForm[graphics, animOpts], 
              CharacterEncoding -> None]];  
          Close[nfilename]; 
ps = Unprotect[Real];  
          Format[x_Real, InputForm] =. ;  
          Protect @@ ps; 
WriteHTMLFile[LiveJarPath, 
            mfilename]];

\(WriteHTMLFile[path_, filename_]\  := \ 
      Module[{nfilename, \ hfilename}, 
hfilename\  = \ 
          StringTake[
              filename, {1, 
                First[First[
                    StringPosition[
                      filename, "."]]]}] <> "html"; 
nfilename = 
          StringJoin[path, hfilename]; 
WriteString[
          nfilename, "<HTML>
"]; 
WriteString[
          nfilename, "<APPLET ARCHIVE='live.jar' CODE='Live.class' WIDTH='650' HEIGHT='650' ALIGN=LEFT VSPACE='50' HSPACE='50'>
"]; 
WriteString[
          nfilename, "<PARAM NAME=BGCOLOR VALUE=#FFFFFF>
"]; 
WriteString[
          nfilename, "<PARAM NAME=MAGNIFICATION VALUE=1.>
"]; 
WriteString[
          nfilename, "<PARAM NAME=INPUT_FILE VALUE='" <> 
            filename <> "'>
"]; 
WriteString[
          nfilename, "</APPLET>
"]; 
WriteString[
          nfilename, "</HTML>"]; 
Close[
          nfilename]; 
WriteLinkButton[path, 
          hfilename]]\)

\(WriteLinkButton[path_, filename_]\  := \ 
      Module[{nb, npath}, 
nb\  = \ 
          InputNotebook[]; 
npath\  = 
          StringReplace[path, {"\" -> "/"}] <> 
            filename; 
NotebookWrite[nb, 
          Cell[TextData[{"View and manipulate the graphic in : ", 
                ButtonBox["browser", 
                  ButtonData -> {URL["file:///" <> npath], None}, 
                  ButtonStyle -> "Hyperlink"], "."}], "Output", CellFrame -> True, FontSize -> 16, 
            FontWeight -> "Bold", 
            Background -> RGBColor[1, 1, 0]]]]\)

Gradient

Unprotect[Gradient];
Gradient[f_] := 
      Module[{temp}, 
temp = 
          Outer[D, {f}, {$xvar, $yvar, $zvar}] // 
            First; 
If[temp[[3]] === 0, Take[temp, 2], 
          temp]]
Protect[Gradient];

\(RotatingDirectionalDerivatives[
        f_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, {x0_, y0_}, n_, 
        opts___?OptionQ] := 
      Module[{thick = 0.02, epsilon = 0.02, u, su, dd, θ, 
          dθ\  = \ 2  π/n, grd, mag}, 
grd = 
          Gradient[f] /. {xvar -> x0, yvar -> y0} // 
            N; 
mag = \@\(grd[\([1]\)]\^2 + grd[\([2]\)]\^2\); 
grd = 
          If[mag > 0, \(1\/mag\) grd, grd]; 
gr3 = 
          ParametricPlot3D[{x0 + t\ grd[\([1]\)], y0 + t\ grd[\([2]\)], 
                epsilon + \((f /. {xvar -> x0 + \ t\ grd[\([1]\)], 
                        yvar -> \ y0 + t\ grd[\([2]\)]})\), Hue[1]}, {t, 
                0, 1}, DisplayFunction -> 
                Identity] /. {List[Hue[aa_], Line[ft_]]\  -> \ 
                List[Hue[aa], Thickness[thick], 
                  Line[ft]]}; 
Table[
          u = {Cos[θ], Sin[θ]}; 
          su = Sequence @@ u; 
dd = grd . u; 
          gr1 = Plot3D[f, \ {xvar, xmin, xmax}, {yvar, ymin, ymax}, 
              Evaluate[opts], ColorFunction -> Hue, 
              PlotPoints -> 25, \ 
              PlotLabel -> 
                StyleForm[\*""(Subscript[f, u])("" <> ToString[x0] <> "," <> 
                    ToString[y0] <> ") = " <> ToString[dd], 
                  FontSize -> 14], 
              DisplayFunction -> Identity]; 
gr2 = 
            ParametricPlot3D[{x0 + t\ Cos[θ], y0 + t\ Sin[θ], 
                  epsilon + \((f /. {xvar -> \ x0 + t\ Cos[θ], 
                          yvar -> \ y0 + t\ Sin[θ]})\)}, {t, 0, 
                  1}, DisplayFunction -> 
                  Identity] /. {List[Line[ft_], rest__]\  -> \ 
                  List[Thickness[thick], Line[ft], 
                    rest]}; 
\(Show[{gr1, gr3, gr2}, 
              DisplayFunction -> $DisplayFunction];\), {θ, 0, 
            2  π - dθ, dθ}];
]\)

Unprotect[Hessian];
Hessian[f_, {x_}] := {D[f, x, x]}
Hessian[
        f_, {x_, y_}] := {{D[f, x, x], D[f, x, y]}, {D[f, y, x], 
          D[f, y, y]}}
Hessian[
        f_, {x_, y_, z_}] := {{D[f, x, x], D[f, x, y], 
          D[f, x, z]}, {D[f, y, x], D[f, y, y], D[f, y, z]}, {D[f, z, x], 
          D[f, z, y], D[f, z, z]}}
Protect[Hessian];

LinearApproximation[f_, {vars__}][
        point__] := (f /. 
            MapThread[Rule[#1, #2] &, {{vars}, {point}}]) + Gradient[
              f, {vars}][point] . 
          First[{{vars} - {point}}]

        
QuadraticApproximation[f_, {vars__}][
        point__] := LinearApproximation[f, {vars}][point] + 
        1/2 (Hessian[f, {vars}] /. 
                MapThread[Rule[#1, #2] &, {{vars}, {point}}]) . 
            First[{{vars} - {point}}] . First[{{vars} - {point}}]

centered functions

\(zoomFunction[{tag_, multiplier_}] := 
      Module[{nb = EvaluationNotebook[], cell, \ 
          rad}, 
NotebookFind[nb, tag, Previous, 
          CellTags]; 

        SelectionMove[nb, All, CellContents]; 
cell\  = \ 
          NotebookRead[nb]; 
\ 
        cell = \ cell /. 
            RowBox[{"Radii", "->", x_}] :> 
              RowBox[{"Radii", "->", 
                  ToBoxes[
                    ToExpression[x]*
                      multiplier]}]\ ; 
cell = \ 
          cell /. RowBox[{"Radii", "->", x_}] :> 
              RowBox[{"Radii", "->", 
                  ToBoxes[
                    ToExpression[x]*
                      multiplier]}]\ ; 
NotebookWrite[nb, 
          cell]; SelectionMove[nb, All, Cell]; 
        SelectionEvaluate[nb, None]
]\)

Options[CenteredContourPlot] = {Zoomable -> False,  
          ColorFunction -> Hue,  ShowGrid -> True,  
          Radii -> 1};

\(CenteredContourPlot[g_, \ {xvar_, x0_}, {yvar_, y0_}, \ 
        opts___?OptionQ] := 
      Module[{zoom, \ contopts, \ gridQ, \ grid, k, \ contONLYopts, \ rad, \ 
          dx, dy, \ gr, nb, cont, rnd, \ 
          semicolonQ}, 
rad = 
          Radii /. Flatten[{opts, 
                Options[CenteredContourPlot]}]; 
zoom = 
          Zoomable /. 
            Flatten[{opts, 
                Options[CenteredContourPlot]}]; 
gridQ = 
          ShowGrid /. 
            Flatten[{opts, 
                Options[CenteredContourPlot]}]; 
If[
          Head[rad] === List, dx = First[rad]; dy = Last[rad], 
          dx = \(dy = 
              rad\)]; 
grid\  = {Epilog -> {RGBColor[
                  0, 0, 0], 
                Table[Line[{{k, y0 - dy}, {k, y0 + dy}}], {k, x0 - dx, \ 
                    x0 + dx, \ dx/5}], \ 
                Table[Line[{{x0 - dx, k}, {x0 + dx, k}}], {k, y0 - dy, \ 
                    y0 + dy, \ dy/5}], PointSize[0.03], 
                Point[{x0, y0}]}}; 
contopts\  = {\ 
            FilterOptions[ContourPlot, opts], 
            ColorFunction -> Hue}; 
contONLYopts = 
          Complement[contopts, 
            Options[CenteredContourPlot]]; 
 (*\(Print["radius = " <> ToString[rad] <> "
zoom = " <> 
                ToString[zoom] <> "
grid = " <> 
                ToString[gridQ] <> "
contopts = " <> 
                ToString[contopts] <> "
contONLYopts = " <> 
                ToString[contONLYopts] <> "
dx = " <> 
                ToString[dx] <> "
dy = " <> 
                ToString[dy]];\)*) 
If[gridQ, \ 
          contopts\  = \ Flatten[{grid, contopts}]]; 
        gr = ContourPlot[
            g, {xvar, x0 - dx, x0 + dx}, {yvar, y0 - dy, y0 + dy}, 
            Evaluate[Apply[Sequence, contopts]]]; 

        If[zoom, {rnd, semicolonQ} = 
            CenteredContourPlotZoom[]; 
If[\(! semicolonQ\), 
\ 
            NotebookWrite[EvaluationNotebook[], 
              Cell[DisplayForm[ToBoxes[gr]], \ "Output", \ 
                CellLabel -> "Out[" <> ToString[$Line] <> "]=", \ 
                CellTags -> "ccp" <> 
                    ToString[
                      rnd] <> "out"]]]]; 
\ \ gr
]\)

\(CenteredContourPlotZoom[] := 
      Module[{nb, cont, rnd, d, \ semiQ}, 
rnd = 
          Random[Integer, {0, 10000}]; 
nb = 
          InputNotebook[]; 
NotebookFind[
          nb, "In[" <> ToString[$Line] <> "]:=", Previous, 
          CellLabel]; 
SelectionMove[nb, After, \ 
          CellContents]; 
NotebookFind[nb, "]", 
          Previous, CellContents]; 
\ 
        SelectionMove[nb, Next, \ Character]; \ 
cont = 
          NotebookRead[nb]; 
SelectionMove[nb, After, \ 
          Character]; 
semiQ\  = 
          If[cont === {}, False, \ True]; 
NotebookFind[
          nb, "In[" <> ToString[$Line] <> "]:=", Previous, 
          CellLabel]; \ 
\ 
        NotebookWrite[nb, 
          NotebookRead[nb] /. Cell[b__] -> 
              Cell[b, CellTags -> "ccp" <> 
                    ToString[rnd]]]; 
SelectionMove[nb, 
          After, Cell]; 
NotebookWrite[nb, 
          Cell[BoxData[
              RowBox[{ButtonBox["zoom out", \ \ ButtonEvaluator -> 
                      Automatic, 
                    ButtonData -> {"ccp" <> ToString[rnd], 2}, \ 
                    ButtonFunction :> \((zoomFunction[#] &)\)], 
                  ButtonBox["zoom in", \ ButtonEvaluator -> Automatic, 
                    ButtonData -> {"ccp" <> ToString[rnd], 1/2}, 
                    ButtonFunction :> \((zoomFunction[#] &)\)]}]], \ "Output", Active -> True, 
            CellTags -> "ccp" <> 
                ToString[
                  rnd] <> "out"]]; 
\ \ SelectionMove[
          nb, Next, Cell]; \ 
        NotebookWrite[nb, 
          NotebookRead[nb] /. Cell[b__] -> 
              Cell[b, CellTags -> "ccp" <> 
                    ToString[rnd] <> "out"]]; \ 
        success = 
          NotebookFind[nb, "Out[" <> ToString[$Line] <> "]=", 
            Previous, CellLabel]; \ 
{rnd, \ 
          semiQ}]\)

Unprotect[ContourPlot];

ContourPlot[
            f_, {xvar_, xmin_,  xmax_}, {yvar_, ymin_, ymax_},  
            Focus ->  {x0_, y0_},  Radii -> {rx_, ry_}, opts___?OptionQ] := 
CenteredContourPlot[
          f, {xvar, x0}, {yvar,  y0}, Radii -> {rx,  ry}, opts]

ContourPlot[
            f_, {xvar_, xmin_,  xmax_}, {yvar_, ymin_, ymax_},  
            Focus ->  {x0_, y0_},  Radii -> rx_, opts___?OptionQ] := 
CenteredContourPlot[
          f, {xvar, x0}, {yvar,  y0}, Radii -> rx, opts]

Protect[
        ContourPlot];

Options[
          CenteredNumericPlot3D] = {Zoomable -> False,   Radii -> 
            1, Steps -> 2};

\(CenteredNumericPlot3D[f_, {xvar_, x0_}, {yvar_, y0_}, 
          opts___?OptionQ]\  := \ 
        Module[{dx, dy, zoom, rad, \ rnd, steps}, 
rad = 
            Radii /. Flatten[{opts, 
                  Options[
                    CenteredNumericPlot3D]}]; 
steps = 
            Steps /. Flatten[{opts, 
                  Options[CenteredNumericPlot3D]}]; 
zoom = 
            Zoomable /. 
              Flatten[{opts, 
                  Options[CenteredNumericPlot3D]}]; 
If[
            Head[rad] === List, dx = First[rad]; dy = Last[rad], 
            dx = \(dy = rad\)]; 
If[\(! zoom\), 
            DisplayForm[\ \ GridBox[
                focusedtableofvalues[
                  f, \ {xvar, \ x0}, {yvar, y0}, \ {dx/steps, dy/steps}], 
                GridFrame -> True, RowLines -> True, 
                RowSpacings -> {2, 0}, ColumnSpacings -> {1, 0.5}, 
                ColumnLines -> True]], 
rnd = 
              Random[Integer, {0, 10000}]; 
nb = 
              InputNotebook[]; 
NotebookFind[
              nb, "In[" <> ToString[$Line] <> "]:=", Previous, 
              CellLabel]; \ 
\ 
            NotebookWrite[nb, 
              NotebookRead[nb] /. Cell[b__] -> 
                  Cell[b, CellTags -> "cnp" <> ToString[rnd]]]; 

            SelectionMove[nb, After, Cell]; 

            NotebookWrite[nb, 
              Cell[BoxData[
                  RowBox[{ButtonBox["zoom out", \ \ ButtonEvaluator -> 
                          Automatic, 
                        ButtonData -> {"cnp" <> ToString[rnd], 2}, \ 
                        ButtonFunction :> \((zoomFunction[#] &)\)], 
                      ButtonBox["zoom in", \ 
                        ButtonEvaluator -> Automatic, 
                        ButtonData -> {"cnp" <> ToString[rnd], 1/2}, 
                        ButtonFunction :> \((zoomFunction[#] &)\)]}]], \ "Output", Active -> True, 
                CellTags -> "cnp" <> 
                    ToString[rnd] <> "out"]]; 
\ \ SelectionMove[nb, 
              After, Cell]; \ 
            NotebookWrite[
              nb, Cell[
                DisplayForm[\ \ GridBox[
                    focusedtableofvalues[
                      f, \ {xvar, \ x0}, {yvar, y0}, \ {dx, dy}], 
                    GridFrame -> True, RowLines -> True, 
                    RowSpacings -> {2, 0}, 
                    ColumnSpacings -> {1, 0.5}, 
                    ColumnLines -> True]], "Output", 
                CellTags -> "cnp" <> 
                    ToString[rnd] <> "out"]];\ \ 
]]\)

CenteredPlot3D // Clear

Options[CenteredPlot3D] = {Zoomable -> False,  
          ColorFunction -> Hue,  Tracings -> True,  
          Radii -> 1};

\(CenteredPlot3D[g_, \ {xvar_, x0_}, {yvar_, y0_}, \ opts___?OptionQ] := 
      Module[{zoom, \ contopts, \ gridQ, \ grid, k, \ contONLYopts, \ rad, \ 
          dx, dy, \ gr, nb, cont, rnd, \ semicolonQ, \ gr1, \ 
          gr2}, 
rad = 
          Radii /. Flatten[{opts, 
                Options[CenteredPlot3D]}]; 
zoom = 
          Zoomable /. 
            Flatten[{opts, 
                Options[CenteredContourPlot]}]; 
gridQ = 
          Tracings /. 
            Flatten[{opts, 
                Options[CenteredContourPlot]}]; 
If[
          Head[rad] === List, dx = First[rad]; dy = Last[rad], 
          dx = \(dy = rad\)]; 
contopts\  = {\ 
            FilterOptions[CenteredPlot3D, opts], 
            FilterOptions[Plot3DExt, opts], 
            Options[CenteredPlot3D]}; 
 (*\(Print["radius = " <> ToString[rad] <> "
zoom = " <> 
                ToString[zoom] <> "
grid = " <> 
                ToString[gridQ] <> "
contopts = " <> 
                ToString[contopts] <> "
contONLYopts = " <> 
                ToString[contONLYopts] <> "
dx = " <> 
                ToString[dx] <> "
dy = " <> 
                ToString[dy]];\)*) 
If[gridQ, \ 
          contopts\  = \ Flatten[{grid, contopts}]]; 
gr1 = 
          Plot3DExt[g, {x, x0 - dx, x0 + dx}, {y, y0 - dy, y0 + dy}, 
            Focus -> {x0, y0}, Evaluate[Apply[Sequence, contopts]], 
            DisplayFunction -> Identity]; 
		gr2\  = \ 
          Graphics3D[{RGBColor[0, 0, 0], PointSize[0.03], \ 
              Point[{x0, y0, g /. {xvar -> x0, yvar -> y0}}]}, 
            DisplayFunction -> Identity]; 
	
        gr = Show[{gr1, gr2}, DisplayFunction -> $DisplayFunction]; 

        If[zoom, {rnd, semicolonQ} = 
            CenteredPlot3DZoom[]; 
If[\(! semicolonQ\), 
\ 
            NotebookWrite[EvaluationNotebook[], 
              Cell[DisplayForm[ToBoxes[gr]], \ "Output", \ 
                CellLabel -> "Out[" <> ToString[$Line] <> "]=", \ 
                CellTags -> "cp" <> 
                    ToString[
                      rnd] <> "out"]]]]; 
\ \ gr
]\)

\(CenteredPlot3DZoom[] := 
      Module[{nb, cont, rnd, d, \ semiQ}, 
rnd = 
          Random[Integer, {0, 10000}]; 
nb = 
          InputNotebook[]; 
NotebookFind[
          nb, "In[" <> ToString[$Line] <> "]:=", Previous, 
          CellLabel]; 
SelectionMove[nb, After, \ 
          CellContents]; 
NotebookFind[nb, "]", 
          Previous, CellContents]; 
\ 
        SelectionMove[nb, Next, \ Character]; \ 
cont = 
          NotebookRead[nb]; 
SelectionMove[nb, After, \ 
          Character]; 
semiQ\  = 
          If[cont === {}, False, \ True]; 
NotebookFind[
          nb, "In[" <> ToString[$Line] <> "]:=", Previous, 
          CellLabel]; \ 
\ 
        NotebookWrite[nb, 
          NotebookRead[nb] /. Cell[b__] -> 
              Cell[b, CellTags -> "cp" <> 
                    ToString[rnd]]]; 
SelectionMove[nb, 
          After, Cell]; 
NotebookWrite[nb, 
          Cell[BoxData[
              RowBox[{ButtonBox["zoom out", \ \ ButtonEvaluator -> 
                      Automatic, 
                    ButtonData -> {"cp" <> ToString[rnd], 2}, \ 
                    ButtonFunction :> \((zoomFunction[#] &)\)], 
                  ButtonBox["zoom in", \ ButtonEvaluator -> Automatic, 
                    ButtonData -> {"cp" <> ToString[rnd], 1/2}, 
                    ButtonFunction :> \((zoomFunction[#] &)\)]}]], \ "Output", Active -> True, 
            CellTags -> "cp" <> 
                ToString[
                  rnd] <> "out"]]; 
\ \ SelectionMove[
          nb, Next, Cell]; \ 
        NotebookWrite[nb, 
          NotebookRead[nb] /. Cell[b__] -> 
              Cell[b, CellTags -> "cp" <> 
                    ToString[rnd] <> "out"]]; \ 
        success = 
          NotebookFind[nb, "Out[" <> ToString[$Line] <> "]=", 
            Previous, CellLabel]; \ 
{rnd, \ 
          semiQ}]\)

NumericPlot3D

FormatEntry[x_]  :=  
      If[x === DirectedInfinity[], x = "∞", 
          If[Head[x] === Integer  ||  Head[x] === Rational, x, 
            PaddedForm[x, {6, 3}]]] // Chop

PlainFormatEntry[x_]  :=  
      If[Head[x] === Integer  ||  Head[x] === Rational, x, N[x, 3]] // 
        Chop

tableofvalues[
          f_, {xvar_, a_, b_, c_:  0.01}, {yvar_, d_, e_, g_:  0.01}] :=Module[{temp,  headings},
\(temp = 
                Table[FormatEntry[f /. {xvar -> i, yvar -> j}] // 
                    ToBoxes, \ {j, \ e, d, \ \(-g\)}, 
	{i, a, b, c}]\);
headings  = 
              Table[ButtonBox[ToBoxes[i], 
                  Background -> RGBColor[1, 1, 0]], {i, a, b, c}];
temp  =  Join[{headings}, temp];
headings  =  
              Join[{ButtonBox["y\x", Background -> RGBColor[1, 0, 1]]}, 
                Table[ButtonBox[ToBoxes[j], 
                    Background -> RGBColor[1, 1, 0]], {j,  e, 
                    d,  -g}]];
temp  =  
              Join[{headings}, Transpose[temp]] // Transpose]

plaintableofvalues[f_, {xvar_, a_, b_, c_:  0.01}, {yvar_, d_, e_, 
            g_:  0.01}] :=Module[{temp,  headings},
\(temp = 
                Table[PlainFormatEntry[f /. {xvar -> i, yvar -> j}], \ {j, \ 
                    e, d, \ \(-g\)}, 
	{i, a, b, c}]\);
headings  = Table[i, {i, a, b, c}];
temp  =  Join[{headings}, temp];
headings  =  
              Join[{"y\x"}, Table[j, {j,  e, d,  -g}]];
temp  =  
              Join[{headings}, Transpose[temp]] // Transpose]

focusedtableofvalues[
          f_, {xvar_, x0_}, {yvar_, y0_}, {dx_, dy_}] :=Module[{temp,  headings},
\(temp = 
                Table[FormatEntry[f /. {xvar -> i, yvar -> j}] // 
                    ToBoxes, \ {j, \ y0 + 2  dy, 
                    y0 - 2  dy, \ \(-dy\)}, 
	{i, x0 - 2  dx, x0 + 2  dx, 
                    dx}]\);
headings  = 
              Table[ButtonBox[ToBoxes[i], 
                  Background -> RGBColor[1, 1, 0]], {i, x0 - 2  dx, 
                  x0 + 2  dx, dx}];
		headings[[3]]  =  
              headings[[3]] /. RGBColor[1, 1, 0] -> RGBColor[0, 1, 1];
		temp[[3, 3]] = 
              ButtonBox[temp[[3, 3]], Background -> RGBColor[0, 1, 1]];
temp  =  Join[{headings}, temp];
headings  =  
              Join[{ButtonBox["y\x", Background -> RGBColor[1, 0, 1]]}, 
                Table[ButtonBox[ToBoxes[j], 
                    Background -> RGBColor[1, 1, 0]], {j, y0 + 2  dy, 
                    y0 - 2  dy,  -dy}]];
		headings[[4]]  =  
              headings[[4]] /. RGBColor[1, 1, 0] -> RGBColor[0, 1, 1];
temp  =  
              Join[{headings}, Transpose[temp]] // Transpose]

NumericPlot3D[
          f_, {xvar_, a_, b_, c_:  1}, {yvar_, d_, e_, g_:  1}] := 
DisplayForm[  GridBox[
          tableofvalues[f,  {xvar,  a,  b,  c}, {yvar,  d,  e,  g}], 
          GridFrame -> True, RowLines -> True, 
          RowSpacings -> {2, 0}, ColumnSpacings -> {1, 0.5}, 
          ColumnLines -> True]]

NumericPlot3D[f_, {xvar_, a_, b_, c_ : 1}, {yvar_, d_, e_, g_ : 1}, 
        Frame -> False] := 
      PlainNumericPlot3D[f, {xvar, a, b, c}, {yvar, d, e, g}]

NumericPlot3D[f_, {xvar_, a_, b_, c_ : 1}, {yvar_, d_, e_, g_ : 1}, 
        Focus -> {x0_, y0_}, Radii -> {rx_, ry_}] := 
      CenteredNumericPlot3D[f, {xvar, x0}, {yvar, y0}, 
        Radii -> {rx, ry}]

NumericPlot3D[f_, {xvar_, a_, b_, c_ : 1}, {yvar_, d_, e_, g_ : 1}, 
        Focus -> {x0_, y0_}, Radii -> rx_] := 
      CenteredNumericPlot3D[f, {xvar, x0}, {yvar, y0}, 
        Radii -> rx]

PlainNumericPlot3D[f_, {xvar_, a_, b_, c_:  1}, {yvar_, d_, e_, g_:  1}] := 
TableForm[
        plaintableofvalues[
          f,  {xvar,  a,  b,  c}, {yvar,  d,  e,  g}]]

Plot3dalt

PlotRangeUsed[
        f_, {xvar_ ,  xmin_, xmax_}, {yvar_, ymin_, ymax_}] :=  
      Module[{gr}, 
gr  =  
          Plot3D[f, {xvar,  xmin, xmax}, {yvar,  ymin, ymax}, 
            DisplayFunction -> 
              Identity]; 
(PlotRange  /.  
              FullOptions[gr])[[3]]
]

TangentPlane[f_, {xvar_, x0_}, {yvar_, y0_}] := z==(f /. {xvar -> x0, yvar -> y0}) + (D[f, x] /. {xvar -> x0}) (x - x0)+(D[f, y] /. {yvar -> y0}) (y - y0) //Expand//Simplify

SetOptions[Plot3D, ColorFunction -> Hue, 
        PlotPoints -> 25,  Mesh -> False, 
        AxesLabel -> {"x", "y", "z"}];

Options[Plot3DExt]  =  Options[Plot3D];
Options[Plot3DExt]  =  
        Join[Options[Plot3DExt], {Focus -> False, 
            TangentLinesAtFocus -> False, 
            PlaneColor -> RGBColor[ .5,  .5,  .5], 
            PlaneEquation ->  False,  
            CrossSectionView -> False,  BothViews -> False, 
            Tracings -> True}];

\(Plot3DExt[g_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, 
        opts___?OptionQ]\  := \ 
      Module[{gr, altgr1, count = 0, \ thick = 0.015, 
          focus = Focus /. Flatten[{opts, Options[Plot3DExt]}], 
          planeColor = 
            PlaneColor /. 
              Flatten[{opts, 
                  Options[Plot3DExt]}], 
viewPoint = 
            ViewPoint /. Flatten[{opts, Options[Plot3DExt]}], 
          planeEquation = 
            PlaneEquation /. Flatten[{opts, Options[Plot3DExt]}], 
          crossSectionView = 
            CrossSectionView /. Flatten[{opts, Options[Plot3DExt]}], 
          bothViews = BothViews /. Flatten[{opts, Options[Plot3DExt]}], 
          plot3dopts = FilterOptions[Plot3D, opts], 
          tangents\  = \ 
            TangentLinesAtFocus /. Flatten[{opts, Options[Plot3DExt]}], 
          paropts = FilterOptions[ParametricPlot3D, opts], 
          graphopts = FilterOptions[Graphics3D, opts], i, 
          rng}, 
gr[\(++count\)] = 
          Plot3D[g, {xvar, xmin, xmax}, {yvar, ymin, ymax}, 
            plot3dopts // Evaluate, 
            DisplayFunction -> Identity]; 
rng\  = 
          PlotRange\  /. \ FullOptions[gr[1]]; 
If[
          VectorQ[focus] && 
            Dimensions[
                focus] == {2}, 
gr[\(++count\)] = 
            Plot3DFocus[g, {xvar, xmin, xmax}, {yvar, ymin, ymax}, focus, \ 
              rng, paropts]; 
If[
            VectorQ[tangents] && Dimensions[tangents] == {2}\  && \ 
              And @@ Map[BooleanQ, tangents], 
If[
              tangents[\([1]\)], 
              gr[\(++count\)] = 
                Plot3Dwxtangent[g, {xvar, xmin, xmax}, {yvar, ymin, ymax}, 
                    focus, \ rng, 
                    paropts] /. {Hue[a_] -> 
                      Sequence[Hue[a], 
                        Thickness[thick]]}]; 
If[
              tangents[\([2]\)], 
              gr[\(++count\)] = 
                Plot3Dwytangent[g, {xvar, xmin, xmax}, {yvar, ymin, ymax}, 
                    focus, \ rng, 
                    paropts] /. {Hue[a_] -> 
                      Sequence[Hue[a], 
                        Thickness[thick]]}]]]; 
If[
          planeEquation =!= False, 
          gr[\(++count\)] = 
            Plot3DWPlane[g, {xvar, xmin, xmax}, {yvar, ymin, ymax}, \ rng, 
              opts]]; 
 (*If[\(crossSectionView === 
              True\)\(,\)]; \ finish\ *) 
Show[
          Table[gr[i], {i, 1, count}], Evaluate[graphopts], 
          ViewPoint -> \ viewPoint, 
          DisplayFunction -> $DisplayFunction, \ 
          AxesLabel -> {"x", "y", "z"}]
]\)

Main3DPlot[g_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, 
        opts___?OptionQ]  := 
      Plot3D[g, {xvar,  xmin, xmax}, {yvar, ymin, ymax}, opts // Evaluate, 
        PlotPoints -> 25,  Mesh -> False,  
        ColorFunction -> Hue,  
        DisplayFunction -> Identity]

Plot3DFocus[
        g_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, {x0_, y0_}, rng_, 
        opts___?OptionQ]  :=  
      Module[{gr2, gr2b, gr3, gr3b,  epsilon, 
          thick = 0.015}, 
epsilon  =  (rng[[3, 
                    2]] - rng[[3, 1]])/
            50; 
gr2  =  
          ParametricPlot3D[{x0, yvar, 
              epsilon + (g /. {xvar -> x0})}, {yvar, ymin, ymax},  
            opts // Evaluate, 
            DisplayFunction -> Identity]; 
gr2b  =  
          ParametricPlot3D[{x0, 
              yvar, -epsilon + (g /. {xvar -> x0})}, {yvar, 
              ymin, ymax},  opts // Evaluate, 
            DisplayFunction -> Identity]; 
gr3  =  
          ParametricPlot3D[{xvar, 
              y0, (g /. {yvar -> y0}) + epsilon}, {xvar, xmin, 
              xmax}, opts // Evaluate, 
            DisplayFunction -> Identity]; 
gr3b  =  
          ParametricPlot3D[{xvar, 
              y0, (g /. {yvar -> y0}) + -epsilon}, {xvar, xmin, 
              xmax}, opts // Evaluate, 
            DisplayFunction -> Identity]; 
Sequence @@ 
          Flatten[{{gr2, gr3, gr2b, 
                  gr3b} /. {List[Line[f_], rest__]  ->  
                    List[Thickness[thick], Line[f], rest]}, 
              Graphics3D[{PointSize[1.5*thick], 
                  Point[{x0, y0, 
                      epsilon + (g /. {xvar -> x0, 
                              yvar -> y0})}], 
                  Point[{x0, 
                      y0, -epsilon + (g /. {xvar -> x0, 
                              yvar -> 
                                y0})}]}]}]
]

Plot3Dwytangent[
        g_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, {x0_, y0_}, rng_, 
        opts___?OptionQ]  :=  
      Module[{epsilon}, 
epsilon  =  (rng[[3, 2]] - 
                rng[[3, 1]])/
            50; 
ParametricPlot3D[{x0, yvar, 
              epsilon + (g /. {xvar -> x0, 
                      yvar -> y0})  +  (D[g /. {xvar -> x0}, 
                        yvar] /. {yvar ->  y0}) (yvar - y0),  
              Hue[1]} // Evaluate, {yvar, ymin, ymax},  opts // Evaluate, 
          DisplayFunction -> Identity]
]

Plot3Dwxtangent[
        g_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, {x0_, y0_}, rng_, 
        opts___?OptionQ]  :=  
      Module[{ epsilon}, 
epsilon  =  (rng[[3, 
                    2]] - rng[[3, 1]])/
            50; 
ParametricPlot3D[{xvar, y0, 
              epsilon + (g /. {xvar -> x0, 
                      yvar -> y0})  +  (D[g /. {yvar -> y0}, 
                        xvar] /. {xvar ->  x0}) (xvar - x0),  
              Hue[1]} // Evaluate, {xvar, xmin, xmax},  opts // Evaluate, 
          DisplayFunction -> Identity]
]

varInPlaneEquationQ[eq_, var_] := 
      Flatten[Position[FullForm[eq], var]] =!= {}

VariablesInPlaneEquation[eq_]  := {varInPlaneEquationQ[eq, x], 
        varInPlaneEquationQ[eq, y], varInPlaneEquationQ[eq, z]}

SolveFor[eq_, var_]  :=  var /. Solve[eq, var] // First

Plot3DWPlane[g_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_},  rng_, 
        opts___?OptionQ]  :=  
      Module[{ plcolor, plane,  vars,  
          paropts = 
            FilterOptions[ParametricPlot3D, 
              opts]}, 
plcolor  = 
          PlaneColor /. 
            Flatten[{opts, Options[Plot3DExt]}]; 
plane  = 
          PlaneEquation /. 
            Flatten[{opts, Options[Plot3DExt]}]; 
vars  = 
          VariablesInPlaneEquation[plane]; 
If[
          vars[[1]], 
If[
            vars[[2]], 
If[
              vars[[3]], 
  
ParametricPlot3D[{xvar, yvar, 
                  SolveFor[plane, z],  plcolor}, {xvar, xmin, xmax}, {yvar, 
                  ymin, ymax}, PlotPoints -> 4,  Evaluate[ paropts], 
                DisplayFunction -> Identity], 
  
ParametricPlot3D[{xvar, SolveFor[plane, yvar], z,  
                  plcolor}, {xvar, xmin, xmax}, {z, rng[[3, 1]],  
                  rng[[3, 2]]}, PlotPoints -> 10,  
                Evaluate[ paropts], 
                DisplayFunction -> Identity]], 
If[
              vars[[3]], 
  
ParametricPlot3D[{xvar, yvar, SolveFor[plane, z],  
                  plcolor}, {xvar, xmin, xmax}, {yvar, ymin, ymax}, 
                PlotPoints -> 10,  Evaluate[ paropts], 
                DisplayFunction -> Identity], 
  
ParametricPlot3D[{SolveFor[
                    plane, xvar], yvar, z,  plcolor}, {yvar, ymin, 
                  ymax}, {z, rng[[3, 1]],  rng[[3, 2]]}, 
                PlotPoints -> 10,  Evaluate[ paropts], 
                DisplayFunction -> Identity]]], 
If[
            vars[[2]], 
If[
              vars[[3]], 
  
ParametricPlot3D[{xvar, yvar, SolveFor[plane, z],  
                  plcolor}, {xvar, xmin, xmax}, {yvar, ymin, ymax}, 
                PlotPoints -> 10,  Evaluate[ paropts], 
                DisplayFunction -> Identity], 
  
ParametricPlot3D[{xvar, 
                  SolveFor[plane, yvar], z,  plcolor}, {xvar, xmin, 
                  xmax}, {z, rng[[3, 1]],  rng[[3, 2]]}, 
                PlotPoints -> 10,  Evaluate[ paropts], 
                DisplayFunction -> Identity]], 
If[
              vars[[3]], 
  
ParametricPlot3D[{xvar, yvar, 
                  SolveFor[plane, z],  plcolor}, {xvar, xmin, xmax}, {yvar, 
                  ymin, ymax}, PlotPoints -> 10,  Evaluate[ paropts], 
                DisplayFunction -> 
                  Identity], 
Print["not a plane in x, y, \
z"]]]]
]

Options[
          Plot3DWxPlane] = {PlaneColor -> 
            RGBColor[ .5,  .5,  .5]};

Approximations3D

Options[
          ShowApproximations3D] = {Method -> "midpoint", 
          Divisions -> {5, 5},  DisplayVolume -> False,  
          Color -> 
            RGBColor[1, 164/255, 215/255]}; 
  

    notListQ[x_] := Head[x] =!= List

ShowApproximations3D[{fn1_, fn2_}, {var1_, a_, b_}, {var2_, c_, d_}, 
        opts___?OptionQ] := 
      Module[{k, rlist, llist, randlist, prt, thex, deltax, mycolors, fnx, 
          xxx, method, divisions, gr1, gr2, volume, grlist, daQ, rc, 
          finopts = {opts}}, 
        method = 
          Method /. Flatten[{opts, 
                Options[
                  ShowApproximations3D]}]; 
  divisions = 
          Divisions /. Flatten[{opts, 
                Options[ShowApproximations3D]}]; 
        daQ = DisplayVolume /. Flatten[{opts, 
                Options[ShowApproximations3D]}]; 
        rc = Color /. Flatten[{opts, 
                Options[ShowApproximations3D]}]; 
        finopts = DeleteCases[finopts, Color -> _]; 
        finopts = DeleteCases[finopts, DisplayVolume -> _]; 
        finopts = DeleteCases[finopts, Divisions -> _]; 
        finopts = DeleteCases[finopts, Method -> _]; 
        mycolors = 
          rc; {volume, 
            grlist} = {RiemannSum[{fn1, fn2}, {var1, a, b}, {var2, c, d}, 
              Divisions -> divisions, Method -> method, 
              DataToo -> False], 
            ApproxCuboids[{fn1, fn2}, {var1, a, b}, {var2, c, d}, 
              Divisions -> divisions, Method -> method, 
              Color -> rc]};  
        gr1 = Show[Graphics3D[grlist], DisplayFunction -> Identity, 
            Lighting -> False];  
gr2 = 
          Plot3D[fn2, {var1, a, b}, {var2, c, d}, ColorFunction -> Hue, 
            DisplayFunction -> Identity]; 
        gr1 = Show[{gr2, gr1, gr2}, DisplayFunction -> $DisplayFunction, 
            Lighting -> False, finopts]; 
If[daQ, 
          volume, gr1]]

ShowApproximations3D[fn_, {var1_, a_, b_}, {var2_, c_, d_}, 
        opts___?OptionQ] := 
      ShowApproximations3D[{0, fn}, {var1, a, b}, {var2, c, d}, 
        opts]

lowerLeftPoint[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := {a, 
        c}
lowerRightPoint[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := {b, 
        c}
upperLeftPoint[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := {a, 
        d}
upperRightPoint[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := {b, 
        d}
midpointPoint[
        f_, {xvar_, a_, b_}, {yvar_, c_, 
          d_}] := {((a + b)/2), ((c + d)/2)}
maxPoint[
        f_, {xvar_, a_, b_}, {yvar_, c_, 
          d_}] := Flatten[
              Table[{f /. {xvar -> j, yvar -> k}, j, k}, {j, a, 
                  b, ((b - a)/(Max[d - c, b - a]*15))}, {k, c, 
                  d, ((d - c)/(Max[d - c, b - a]*15))}], 1] // Sort // 
          Last // Rest
minPoint[
        f_, {xvar_, a_, b_}, {yvar_, c_, 
          d_}] := Flatten[
              Table[{f /. {xvar -> j, yvar -> k}, j, k}, {j, a, 
                  b, ((b - a)/(Max[d - c, b - a]*15))}, {k, c, 
                  d, ((d - c)/(Max[d - c, b - a]*15))}], 1] // Sort // 
          First // Rest
closestCornerPoint[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := 
      Which[a >= 0 && c >= 0, {a, c}, 
a >= 0 && 
          d <= 0, {a, d}, 
a >= 0 && 
          Sign[c] =!= Sign[d], {a, 
          If[d > Abs[c], c, d]}, 
b <= 0 && c >= 0, {b, 
          c}, 
b <= 0 && d <= 0, {b, 
          d}, 
b <= 0 && Sign[c] =!= Sign[d], {b, 
          If[d > Abs[c], c, d]}, 
c >= 0 && 
          Sign[a] =!= Sign[b], {If[b > Abs[a], a, b], 
          c}, 
d <= 0 && 
          Sign[a] =!= Sign[b], {If[b > Abs[a], a, b], 
          d}, 
Sign[c] =!= Sign[d] && 
          Sign[a] =!= Sign[b], {If[b > Abs[a], a, b], 
          If[d > Abs[c], c, d]}
]

lowerLeftValue[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := 
      f /. {xvar -> a, yvar -> c}
lowerRightValue[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := 
      f /. {xvar -> b, yvar -> c}
upperLeftValue[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := 
      f /. {xvar -> a, yvar -> d}
upperRightValue[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := 
      f /. {xvar -> b, yvar -> d}
midpointValue[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := 
      f /. {xvar -> ((a + b)/2), 
          yvar -> ((c + d)/2)}
minValue[
        f_, {xvar_, a_, b_}, {yvar_, c_, 
          d_}] := 
ApproxValues[
            f, {xvar, a, b}, {yvar, c, d}, Divisions -> {15, 15}, 
            Method -> "midpoint"] // Flatten // 
        Min
maxValue[
        f_, {xvar_, a_, b_}, {yvar_, c_, 
          d_}] := ApproxValues[f, {xvar, a, b}, {yvar, c, d}, 
            Divisions -> {15, 15}, Method -> "midpoint"] // 
          Flatten // Max
closestCornerValue[f_, {xvar_, a_, b_}, {yvar_, c_, d_}] := 
      Which[a >= 0 && c >= 0, 
        f /. {xvar -> a, yvar -> c}, 
a >= 0 && 
          d <= 0, f /. {xvar -> a, 
            yvar -> d}, 
a >= 0 && 
          Sign[c] =!= Sign[d], 
        f /. {xvar -> a, 
            yvar -> If[d > Abs[c], c, d]}, 
b <= 0 && 
          c >= 0, f /. {xvar -> b, 
            yvar -> c}, 
b <= 0 && d <= 0, 
        f /. {xvar -> b, yvar -> d}, 
b <= 0 && 
          Sign[c] =!= Sign[d], 
        f /. {xvar -> b, 
            yvar -> If[d > Abs[c], c, d]}, 
c >= 0 && 
          Sign[a] =!= Sign[b], 
        f /. {xvar -> If[b > Abs[a], a, b], 
            yvar -> c}, 
d <= 0 && 
          Sign[a] =!= Sign[b], 
        f /. {xvar -> If[b > Abs[a], a, b], 
            yvar -> d}, 
Sign[c] =!= Sign[d] && 
          Sign[a] =!= Sign[b], 
        f /. {xvar -> If[b > Abs[a], a, b], 
            yvar -> 
              If[d > Abs[c], c, d]}
]

ApproxPoints[f_, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_] := 
      With[{dx = ((b - a)/m), dy = ((d - c)/n)}, 
        Table[ToExpression[method <> "Point"][
            f, {xvar, a + j dx, a + (j + 1) dx}, {yvar, c + k dy, 
              c + (k + 1) dy}], {j, 0, m - 1}, {k, 0, 
            n - 1}]]

ApproxValues[f_, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_] := 
      With[{dx = ((b - a)/m), dy = ((d - c)/n)}, 
        Table[ToExpression[method <> "Value"][
            f, {xvar, a + j dx, a + (j + 1) dx}, {yvar, c + k dy, 
              c + (k + 1) dy}], {j, 0, m - 1}, {k, 0, 
            n - 1}]]

NApproxValues[f_, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_] := 
      With[{dx = ((b - a)/m), dy = ((d - c)/n)}, 
        Table[ToExpression[method <> "Value"][
            f, {xvar, a*(1.0) + j dx, a + (j + 1) dx}, {yvar, 
              c*(1.0) + k dy, c + (k + 1) dy}], {j, 0, m - 1}, {k, 
            0, n - 1}]]

ApproxCuboids[{f_, g_}, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_, opts___] := 
      Module[{dx = ((b - a)/m), dy = ((d - c)/n), lowvals, uppervals, 
          lowers, uppers, cuboids, color}, 
color = 
          Color /. Flatten[{opts, 
                Options[ShowApproximations3D]}]; 
lowers = 
          Flatten[Table[{a + j dx, c + k dy}, {j, 0, m - 1}, {k, 0, 
                n - 1}], 1]; 
lowvals = 
          If[method === "min", 
Flatten[
              NApproxValues[f, {xvar, a, b}, {yvar, c, d}, 
                Divisions -> {m, n}, Method -> "max"], 1], 
            If[method === "max", 
              Flatten[NApproxValues[f, {xvar, a, b}, {yvar, c, d}, 
                  Divisions -> {m, n}, Method -> "min"], 1], 
              Flatten[NApproxValues[f, {xvar, a, b}, {yvar, c, d}, 
                  Divisions -> {m, n}, Method -> method], 
                1]]]; 
lowers = 
          Join[Transpose[lowers], {lowvals}] // 
            Transpose; 
uppers = 
          Flatten[Table[{a + j dx, c + k dy}, {j, 1, m}, {k, 1, n}], 
            1]; 
uppervals = 
          Flatten[NApproxValues[g, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, Method -> method], 
            1]; 
uppers = 
          Join[Transpose[uppers], {uppervals}] // 
            Transpose; 
cuboids = 
          Transpose[{lowers, uppers}]; 
        Prepend[(Cuboid[First[#1], Last[#1]] &) /@ cuboids, 
          color]]

Options[RiemannSum] = {DataToo -> False,  
          Method -> "midpoint"};

RiemannSum[{lower_, upper_}, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Method -> method_, Divisions -> {m_, n_}, opts___?OptionQ] := 
      RiemannSum[{lower, upper}, {xvar, a, b}, {yvar, c, d}, 
        Divisions -> {m, n}, Method -> method, opts]

\(\(\(RiemannSum[{lower_, upper_}, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_, opts___?OptionQ] := 
      Module[{bottomvals, topvals, xyPoints, data, datatoo, sum, meth, 
          dx = \(b - a\)\/m, dy = \(d - c\)\/n}, 
datatoo = 
          DataToo /. 
            Flatten[{opts, Options[RiemannSum]}]; 
 (*\ 
          Return[{method, opts, Options[RiemannSum], \ 
              Head[method]}]\ *) ; 
Switch[
          method, "min", 
          bottomvals = 
            ApproxValues[lower, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, Method -> "max"]; 

          topvals = 
            ApproxValues[upper, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, 
              Method -> "min"], 
"max", 
          bottomvals = 
            ApproxValues[lower, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, Method -> "min"]; 

          topvals = 
            ApproxValues[upper, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, 
              Method -> "max"], 
_, 
          bottomvals = 
            ApproxValues[lower, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, Method -> method]; 

          topvals = 
            ApproxValues[upper, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, 
              Method -> method]]; 
data\  = \ 
          topvals - bottomvals; 
sum = 
          Plus @@ \((Flatten[data]*dx*dy)\); 
xyPoints = 
          Table[{j, k}, {j, a, b - dx, dx}, {k, c, d - dy, 
              dy}]; 
data = 
          Flatten[xyPoints, 1] // 
            Transpose; 
data\  = {First[data], Last[data], 
              Flatten[bottomvals], Flatten[topvals]} // 
            Transpose; 
If[datatoo, {sum, data}, \ 
          sum]
]\)\(
\)
    \)\)
RiemannSum[f_, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_, opts___?OptionQ] := 
      RiemannSum[{0, f}, {xvar, a, b}, {yvar, c, d}, 
        Divisions -> {m, n},  Method -> method, 
        opts]

    
RiemannSum[f_, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, opts___?OptionQ] := 
      RiemannSum[{0, f}, {xvar, a, b}, {yvar, c, d}, 
        Divisions -> {m, n},  Method -> midpoint, 
        opts]

    
RiemannSum[{lower_, upper_}, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, opts___?OptionQ] := 
      RiemannSum[{lower, upper}, {xvar, a, b}, {yvar, c, d}, 
        Divisions -> {m, n},  Method -> midpoint, 
        opts]

\(\(\(NRiemannSum[{lower_, upper_}, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_, opts___?OptionQ] := 
      Module[{bottomvals, topvals, xyPoints, data, datatoo, sum, meth, 
          dx = \(b - a\)\/\(m*1.0\), 
          dy = \(d - c\)\/\(n*1.0\)}, 
datatoo = 
          DataToo /. 
            Flatten[{opts, Options[RiemannSum]}]; 
 (*\ 
          Return[{method, opts, Options[RiemannSum], \ 
              Head[method]}]\ *) ; 
Switch[
          method, "min", 
          bottomvals = 
            NApproxValues[lower, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, Method -> "max"]; 

          topvals = 
            NApproxValues[upper, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, 
              Method -> "min"], 
"max", 
          bottomvals = 
            NApproxValues[lower, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, Method -> "min"]; 

          topvals = 
            NApproxValues[upper, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, 
              Method -> "max"], 
_, 
          bottomvals = 
            NApproxValues[lower, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, Method -> method]; 

          topvals = 
            NApproxValues[upper, {xvar, a, b}, {yvar, c, d}, 
              Divisions -> {m, n}, 
              Method -> method]]; 
data\  = \ 
          topvals - bottomvals; 
sum = 
          Plus @@ \((Flatten[data]*dx*dy)\); 
xyPoints = 
          Table[{j, k}, {j, a, b - dx, dx}, {k, c, d - dy, 
              dy}]; 
data = 
          Flatten[xyPoints, 1] // 
            Transpose; 
data\  = {First[data], Last[data], 
              Flatten[bottomvals], Flatten[topvals]} // 
            Transpose; 
If[datatoo, {sum, data}, \ 
          sum]
]\)\(
\)
    \)\)
NRiemannSum[f_, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, Method -> method_, opts___?OptionQ] := 
      NRiemannSum[{0, f}, {xvar, a, b}, {yvar, c, d}, 
        Divisions -> {m, n},  Method -> method, 
        opts]

    
NRiemannSum[f_, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, opts___?OptionQ] := 
      NRiemannSum[{0, f}, {xvar, a, b}, {yvar, c, d}, 
        Divisions -> {m, n},  Method -> midpoint, 
        opts]

    
NRiemannSum[{lower_, upper_}, {xvar_, a_, b_}, {yvar_, c_, d_}, 
        Divisions -> {m_, n_}, opts___?OptionQ] := 
      NRiemannSum[{lower, upper}, {xvar, a, b}, {yvar, c, d}, 
        Divisions -> {m, n},  Method -> midpoint, 
        opts]

ParametricPlot

TraceParametricPlot[{ff_, gg_}, {t_, t0_, t1_}, dt_Integer:  8, 
        opts___] := 
      TraceParametricPlot[{f,  g},  {t,  t0,  t1,  dt},  
        opts]

\(\(\(TraceParametricPlot[{ff_, gg_}, {t_, t0_, t1_, dt_Integer:  8}, 
        opts___] := 
\ \ Module[{ndt, arro, arrsz, rng, xx, tt, f, g}, 
        ndt = N[\((t1 - t0)\)/dt]; 
\ \ \ \ f = Function[t, ff]; \ 
        g = Function[t, 
            gg]; 
\ \ \ \ \ arrsz = 
\ \ \ \ \ \ N[\((Max[
                    Table[Abs[f[tt]], {tt, t0, t1, 
                        N[\((t1 - t0)\)/
                            25]}]] + 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ Max[
                    Table[Abs[g[tt]], {tt, t0, t1, N[\((t1 - t0)\)/25]}]])\)/
              2.0]; 
\ \ \ \ \ arro[{a_, b_}, {c_, d_}, 
            eps_] := 
\ \ \ \ \ \ Polygon[{{a, b} + 
                3*eps*{c - a, d - b}/Sqrt[\((b - d)\)^2 + \((c - a)\)^2], {a, 
                  b} + 
\ \ \ \ \ \ \ \ \ \ \ \ eps/
                    Sqrt[\((b - d)\)^2 + \((c - a)\)^2]*{\((b - d)\), \((c - 
                        a)\)}, {a, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ b} + \((\(-eps\))\)/Sqrt[\((b - d)\)^2 + \((c - a)\)^2]*{\((b - d)\), \((c - 
                        a)\)}}]; 
\ \ \ \ \ rng = 
\ \ \ \ \ \ FullOptions[
\ \ \ \ \ \ \ \ ParametricPlot[{f[xx], g[xx]}, {xx, t0, t1}, opts, 
              PlotRange -> 
                All, 
\ \ \ \ \ \ \ \ \ \ DisplayFunction -> Identity], 
            PlotRange]; 
\ \ \ \ \ Do[
\ \ \ \ \ \ ParametricPlot[{f[xx], 
              g[xx]}, {xx, t0, tt}, 
\ \ \ \ \ \ \ \ PlotRange -> 
              rng + {{0, 
                    0}, {0,  .06*
                      rng[\([2, 2]\)]}}, 
\ \ \ \ \ \ \ \ PlotStyle -> 
              Hue[\((tt + 1)\)/\((dt - 
                      1)\)], 
\ \ \ \ \ \ \ \ Epilog -> {Blue, 
\ \ \ \ \ \ \ \ \ \ \ \ arro[{f[tt], g[tt]}, {f[tt +  .05*arrsz], 
                    g[\((tt +  .05*
                            arrsz)\)]},  .035*
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ arrsz]}, 
\ \ \ \ \ \ \ \ PlotLabel
\ \ \ \ \ \ \ \ \ \  -> 
              StringJoin["t = ", 
\ \ \ \ \ \ \ \ \ \ \ \ ToString[
                  N[If[Abs[tt +  .001] <  .0001, 0, tt +  .001], 
                    3]]], 
\ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction, opts], {tt, t0 -  .001, t1, 
            ndt}]]\)\(\ \ \ \ \ \)\)\)

Options[
          RainbowParametricPlot]  = {IncludeDots -> 
            False};

RainbowParametricPlot[{ff_, gg_}, {var_, min_, max_}, n_Integer : 8, 
        opts___?OptionQ] := 
      Module[{ln, gr, rt, i, pt, j, rng, dx, st, grs, f, g, dotsQ}, 
        dotsQ = IncludeDots /. Flatten[{opts, 
                Options[RainbowParametricPlot]}]; 
If[
          dotsQ, RainbowParametricPlotwDots[{ff, gg}, {var, min, max}, 
            n], 
dx = ((max - min)/n); 
f = Function[var, ff]; g = Function[var, gg]; 
          Do[gr[j] = 
              ParametricPlot[{f[var], g[var]}, {var, min + (j - 1) dx, 
                  min + (j - 1) dx + dx}, 
                DisplayFunction -> Identity, 
                PlotStyle -> Hue[(j/n)], opts], {j, 1, n}]; 
          grs = Table[gr[i], {i, 1, n}]; 
          Show[grs, DisplayFunction -> $DisplayFunction]]]

\(RainbowParametricPlotwDots[{ff_, gg_}, {var_, min_, max_}, 
        n_Integer:  8, 
        opts___] := 
\ \ Module[{ln, gr, rt, i, pt, rng, dx, st, grs, f, g, 
          pts = {}}, 
        dx = \((max - min)\)/n; 
\ \ \ \ \ f = 
          Function[var, ff]; 
\ \ \ \ \ g = 
          Function[var, gg]; 
\ \ \ \ \ Do[
          st = min + \((i - 
                    1)\)\ dx; 
\ \ \ \ \ \ \ AppendTo[
\ \ \ \ \ \ \ \ pts, {PointSize[0.01 + i\ 0.0025], Point[{f[st], g[st]}]}]; 
\ \ \ \ \ \ \ gr[
              i] = 
\ \ \ \ \ \ \ \ ParametricPlot[{f[var], g[var]}, {var, 
                st, st + dx}, 
\ \ \ \ \ \ \ \ \ \ DisplayFunction -> 
                Identity, PlotStyle -> Hue[i/n], opts], {i, 
            1, 
\ \ \ \ \ \ \ \ n}]; \ 
        grs = Table[gr[i], {i, 1, n}]; 
\ \ \ \ \ Show[{grs, 
            Graphics[
              pts]}, 
\ \ \ \ \ \ DisplayFunction -> $DisplayFunction]]\)

LineIntegrals

Distance[pt1_List,  
        pt2_List]  :=  (Sqrt[(pt2 - pt1) . (pt2 - pt1)])

\(LineIntegralIllustration[{P_, Q_}, {ff_, gg_}, {var_, min_, max_}, 
        n_Integer:  8, 
        opts___] := 
\ \ Module[{ln, gr, rt, i, pt, rng, dt, st, grs, f, g, 
          pts, scale, vfarrows, disarrows, widest}, 
        dt = \((max - min)\)/n; 
\ \ \ \ \ f = 
          Function[var, ff]; 
\ \ \ \ \ g = 
          Function[var, gg]; 
pts\  = \ 
          Table[Point[{f[k], g[k]}], {k, min, \ max\ , \ 
              dt}]; 
disarrows = 
          Table[Arrow[First[pts[\([k]\)]], First[pts[\([k + 1]\)]]], {k, 1, 
              Length[pts] - 1}]; 
grs = 
          ParametricPlot[{f[var], g[var]}, {var, min, 
              max}, 
\ \ \ \ \ \ \ \ \ \ opts, \ 
            DisplayFunction -> Identity]; \ 
widest = 
          PlotRange\  /. \ FullOptions[grs]; 
widest\  = 
          Map[\((Last[#] - First[#])\) &, widest] // 
            Max; 
vfarrows = 
          Map[Arrow[First[#], 
                  First[#] + {P /. {x -> #[\([1, 1]\)], 
                          y -> #[\([1, 2]\)]}, 
                      Q /. {x -> #[\([1, 1]\)], 
                          y -> #[\([1, 2]\)]}}] &, pts] // 
            Drop[#, \(-1\)] &; 
\ \ \ \ scale = 
          1/5*widest/\((Map[Apply[Distance, #] &, N[vfarrows]] // 
                  Max)\); 
vfarrows = 
          vfarrows /. 
            Arrow[first_, second_] -> 
              Arrow[first, 
                first + scale*\((second - first)\)]; 
\ \ \ \ \ Show[{grs, 
            Graphics[{PointSize[0.02], pts, RGBColor[1, 0, 0], disarrows, 
                RGBColor[0, 0, 1], 
                vfarrows}]}, 
\ \ \ \ \ \ DisplayFunction -> $DisplayFunction]]\)