Find where f(x,y)=(x^2)+x y-2(y^3)-13x+2  (Subscript[f, x]) and (Subscript[f, y]) are both simultaneously 0.

1.	Define the function:

func=(x^2)+x y-2(y^3)-13x+2;

2.	Create a plot of the surface:

Plot3D[func,{x,-10,10},{y,-10,10},BoxRatios->1,AxesLabel->{"x","y","z"}]

3.	Solve  (Subscript[f, x])=(Subscript[f, y])=0 for x and y:

sol=Solve[{D[func,x]==0,D[func,y]==0},{x,y}]

4.	Find the z coordinate of each point:

zVal=func/.sol

5.	Combine the (x,y) points along with z-values to make points in (ℝ^3):

pts=MapThread[Append,{{x,y}/.sol,zVal}]

6.	Plot the surface and points:

Show[{
Plot3D[func,{x,-10,10},{y,-10,10}],
Graphics3D[{Red, PointSize[.05],Point[pts]}]
},AxesLabel->{"x","y","z"},BoxRatios->{1,1,2}]

Instead of using Point we can use Ellipsoid to make better looking points.

Show[{
Plot3D[func,{x,5,9},{y,-2,2}],
Graphics3D[{Red,Ellipsoid[#,{.1,.1,.2}]&/@pts}]
},AxesLabel->{"x","y","z"},BoxRatios->{1,1,2}]