$$f(x,y)=\frac{2}{x^2+(y-1)^2+1}-\frac{1}{x^2+(y+1)^2+1}$$ will do.
Indeed, this function $f$ has three critical points: a saddle point $(0,\approx-8.54)$, a point of a local minimum at $(0,\approx-1.14)$, and a point of a local maximum at $(0,\approx1.04)$.
Below are the calculations, in Mathematica:
To have noncollinear critical points (again, a saddle, a local minimum, and a local maximum), one can use $$f(x,y)=\frac{2}{x^2/2+x/10+(y-1)^2+1}-\frac{1}{x^2+(y+1)^2+1}$$ instead. Then, with a suitable affine transformation $A$, one can have any given three noncollinear points on the plane to be the three critical points (again, a saddle, a local minimum, and a local maximum) of the function $z\mapsto f(Az)$.