Lagrange Multipliers in 2D

Set up function and solve for the constrained critical points

Clear[f,g,fx,fy,gx,gy,x,y,L,eqn,pts,gradf,gradg]
f[x_,y_]:= (y^2)-4 x* y + 4 (x^2)
g[x_,y_]:= (x^2)+  (y^2) - 1

fx[x_,y_]=(∂_x) f[x,y]
fy[x_,y_]=(∂_y) f[x,y]
gradf={ D[f[x,y],x], D[f[x,y],y]}

gx[x_,y_]=(∂_x) g[x,y]
gy[x_,y_]=(∂_y)g[x,y]
gradg={ D[g[x,y],x], D[g[x,y],y]}

eqns={fx[x,y] == λ* gx[x,y],fy[x,y] ==  λ* gy[x,y],g[x,y]==0}
pts=Solve[eqns,{x,y, λ}]
Print["The method requires that the gradient of g is nonzero at a critical point"]
{gx[x,y], gy[x,y]}/.pts

f[x,y]/.pts
N[f[x,y]/.pts]

graphs

Clear[xlo,ylo,xhi,yhi,zlo,zhi,numCon,t1,t2];
xlo=ylo=-2; xhi=yhi=2;
numCon=25;

Clear[fPlot3D,fPlotGradient,fPlotLevel,
gPlot3Da,gPlot3Db,gPlot3DCurvea,gPlotGradient,gPlot3DCurveb,gPlotImplicit,
cpPlot2D,cpPlot3D];

fPlot3D=Plot3D[f[x,y],{x,xlo,xhi},{y,ylo,yhi},AxesLabel->{" x "," y "," z "}]

zlo=-1;zhi=20;

Solve[g[x,y]==0,y]

gPlot3Da=ParametricPlot3D[{x,-(Sqrt[1-(x^2)]),z},{x,xlo,xhi},{z,zlo,zhi}];
gPlot3Db=ParametricPlot3D[{x,(Sqrt[1-(x^2)]),z},{x,xlo,xhi},{z,zlo,zhi}];

Show[{fPlot3D,gPlot3Da,gPlot3Db},AxesLabel->{" x "," y "," z "}]

gPlot3DCurvea=ParametricPlot3D[{x,(Sqrt[1-(x^2)]),f[x,(Sqrt[1-(x^2)])]},{x,xlo,xhi},{z,zlo,zhi},PlotStyle->{Red,Thickness[.15]}];
gPlot3DCurveb=ParametricPlot3D[{x,-(Sqrt[1-(x^2)]),f[x,-(Sqrt[1-(x^2)])]},{x,xlo,xhi},{z,zlo,zhi},PlotStyle->{Red,Thickness[.15]}];
Show[{fPlot3D,gPlot3DCurvea,gPlot3DCurveb},AxesLabel->{" x "," y "," z "}]

Plot3D[f[x,y],{x,xlo,xhi},{y,ylo,yhi},AxesLabel->{" x "," y "," z "},RegionFunction->Function[{x,y,z},-0.1<g[x,y]<0.1]]

fPlotLevel=ContourPlot[f[x,y],{x,xlo,xhi},{y,ylo,yhi},Contours->numCon,ContourShading->False];
gPlotImplicit=ContourPlot[g[x,y]==0,{x,xlo,xhi},{y,ylo,yhi},Axes->True];
Show[{fPlotLevel,gPlotImplicit},Axes->True]

t1=Table[{x,y,f[x,y]}/.pts[[index]],{index,1,Length[pts]}]
cpPlot3D=ListPointPlot3D[t1,PlotStyle->PointSize[0.03]];
Show[{fPlot3D,cpPlot3D},AxesLabel->{" x "," y "," z "}]

t2=Table[Point[{x,y}/.pts[[index]]],{index,1,Length[pts]}];
cpPlot2D = Show[Graphics[{PointSize[0.02`],t2}],PlotRange->All,Axes->True];
Show[{{fPlotLevel,gPlotImplicit}},Graphics[{PointSize[0.02`],t2}],PlotRange->All,Axes->True]

<<"VectorFieldPlots`"

CoolColor[ z_ ] :=RGBColor[z,1-z,1];
fPlotGradient=VectorFieldPlots`VectorFieldPlot[gradf,{x,xlo,xhi},{y,ylo,yhi}, ColorFunction -> CoolColor]
gPlotGradient=VectorFieldPlots`VectorFieldPlot[gradg,{x,xlo,xhi},{y,ylo,yhi}]

the critical points are on the constraint curve where the gradients of g and of f are parallel.

Show[{fPlotGradient,gPlotGradient,cpPlot2D,gPlotImplicit},PlotRange->All,Axes->True]