Introducimos la función

> g:=(x,y)->y/(1+x^2+y^2);

Solve, resuelve directamente el sistema de las dos derivadas parciales primeras nulas

> solve({diff(g(x,y),x)=0,diff(g(x,y),y)=0},{x,y});

Para clasificar los extremos, construimos la matriz Hessiana

> with(linalg):

> hessian(g(x,y),[x,y]):

> H:=(x,y)->hessian(g(x,y),[x,y]);

Evaluamos la matriz H en los puntos críticos y calculamos los determinantes

> H1:=subs(x=0,y=1,H(x,y));

> det(H1);

> H2:=subs(x=0,y=-1,H(x,y));

> det(H2);