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

a=Plot3D[x Exp[-((x^2)+(y^2))],{x,-2,2},{y,-2,2}]

gradf[x_,y_]=Grad[x Exp[-((x^2)+(y^2))],{x,y}]

Together[%]

Solve[{-1+2 (x^2),x y}=={0,0},{x,y}]

{f[((-1)/(Sqrt[2])),0],f[(1/(Sqrt[2])),0]}

N[%]

b=ParametricPlot3D[{x,0,x Exp[-(x^2)]},{x,-2,2},PlotStyle->{Thick,Blue}];

Show[a,b]

c=ContourPlot[x Exp[-((x^2)+(y^2))],{x,-2,2},{y,-2,2},ContourShading->False,ContourLabels->True];

d=StreamPlot[gradf[x,y],{x,-2,2},{y,-2,2}];

f=Show[c,d,Graphics[{{PointSize[Large],Blue,Thick,Point[{(1/(Sqrt[2])),0}]},{PointSize[Large],Red,Thick,Point[{((-1)/(Sqrt[2])),0}]}}]]

lapf[x_,y_]=Div[gradf[x,y],{x,y}]

e=ContourPlot[ArcSinh[lapf[x,y]],{x,-1.5,1.5},{y,-1.5,1.5},Contours->30,ContourStyle->{LightYellow,Opacity[0.01]}];

Show[e,f]