Del[F_]:={D[F[x,y,z],x],D[F[x,y,z],y],D[F[x,y,z],z]}
Del4[F_]:={D[F[x,y,z,λ],x],D[F[x,y,z,λ],y],D[F[x,y,z,λ],z],D[F[x,y,z,λ],λ]}
Lag[F_,G_]:=Module[{Λ,SOL},
Λ[x_,y_,z_,λ_]=F[x,y,z]+λ G[x,y,z];
SOL={λ,x,y,z,F[x,y,z]}/.Solve[Del4[Λ]=={0,0,0,0}];
TableForm[FullSimplify[SOL],TableHeadings->{None,{"λ","x","y","z","Val"}},TableSpacing->{5,2}]
]
NLag[F_,G_]:=Module[{Λ,SOL},
Λ[x_,y_,z_,λ_]=F[x,y,z]+λ G[x,y,z];
SOL={λ,x,y,z,F[x,y,z]}/.Solve[Del4[Λ]=={0,0,0,0}];
TableForm[N[SOL],TableHeadings->{None,{"λ","x","y","z","Val"}},TableSpacing->{5,2}]
]

F[x_,y_,z_]=9(x^2)+6(y^2)-4(z^2);
G[x_,y_,z_]=(x^2)+(y^2)+6(z^2)-4;
Lag[F,G]

vol[x_,y_,z_]=(2x)(2y)(2z);
el[x_,y_,z_]=18(x^2)+43(y^2)+21(z^2)-12;
Lag[vol,el]

Dist[x_,y_,z_]=(Sqrt[(x^2)+(y^2)+(z^2)]);
S[x_,y_,z_]=x y+2x z-5(Sqrt[5]);
Lag[Dist,S]

Dist[x_,y_,z_]=(Sqrt[(((x-2))^2)+(((y-4))^2)+(z^2)]);
S[x_,y_,z_]=10x+x y+2x (z^2)-10;
ContourPlot3D[S[x,y,z]==0,{x,-2,6},{y,0,8},{z,-4,4}];
Graphics3D[{PointSize[Large],Red,Point[{2,4,0}]}];
Show[%,%%]
NLag[Dist,S]

Dist[x_,y_,z_]=(Sqrt[(((x-2))^2)+(((y-4))^2)+(z^2)]);
S[x_,y_,z_]=-10x+x y+2x (z^2)-10;
NLag[Dist,S]
Dist2D[x_,y_]=Dist[x,y,z]/.Solve[S[x,y,z]==0,z][[1]];
{X,Y}={x,y}/.Solve[D[Dist2D[x,y],x]==0,x][[1]]⋃Solve[D[Dist2D[x,y],y]==0,y][[1]];
Z=z/.Solve[S[X,Y,z]==0,z][[1]];
{X,Y,Z,Dist2D[X,Y]}//N
Z=z/.Solve[S[X,Y,z]==0,z][[2]];
{X,Y,Z,Dist2D[X,Y]}//N

ContourPlot3D[S[x,y,z]==0,{x,-20,20},{y,0,20},{z,-7,8}];
Graphics3D[{PointSize[Large],Red,Point[{2,4,0}],PointSize[Large],Green,Point[{2.42509,4.25,-2.22188}],Point[{2.42509,4.25,2.22188}],Point[{5.0833,11.9672,0}],Point[{-1.47905,3.23892,0}]}];
Show[%,%%]
Plot3D[Dist2D[x,y],{x,0.05,20},{y,0,20},PlotPoints->30];
Plot3D[Dist2D[x,y],{x,-20,-0.05},{y,0,20},PlotPoints->30];
Graphics3D[{PointSize[Large],Red,Point[{X,Y,Dist2D[X,Y]}]}];
Show[%,%%,%%%,PlotRange->All]