> e:=exp(1);tay:=taylor(e^x,x=0,3);poly:=convert(tay,'polynom');
> f:=e^(-x^2-(y+1)^2);
> readlib(mtaylor);
> T:=mtaylor(f,[x=0,y=-1],3);
> L:=2;a:=plot3d(f,x=0-L..0+L,y=-1-L..-1+L,style=patchcontour,shading=z):b:=plot3d(T,x=0-L..0+L,y=-1-L..-1+L,style=patchnogrid):with(plots):display(a,b);
> f:=x^2+y^3-3*x*y;#has 2
> f:=x^3-3*x*y+y^3;#has 2 real, 2 complex
> f:=x*y*(3-x-y);#has four
> f:=4*x^2+y^2-4*x+2*y;#has 1
> f:=x*sin(y);#has infinity many but ...
> diff(f,x);diff(f,y);
> critical:=solve({diff(f,x)=0,diff(f,y)=0});
> fxx:=diff(f,x,x);fyy:=diff(f,y,y);fxy:=diff(f,x,y);
> discriminant:=fxx*fyy-fxy*fxy;
> critical[1];#only works with more than one answer
> critical[1];subs(critical[1],[discriminant,fxx]);
> critical;dd:=subs(critical,[discriminant,fxx]);#only works with exactly one answer
> evalf(dd);
> critical[2];subs(critical[2],[discriminant,fxx]);
> critical[3];subs(critical[3],[discriminant,fxx]);
> critical[4];subs(critical[4],[discriminant,fxx]);
> for i from 1 to nops([critical]) do critical[i];subs(critical[i],[discriminant,fxx]);od;#only works when there is more than one answer
> nops([critical]);
> plot3d(f,x=-4..4,y=-4..4);
>