> restart;
> f := exp(-(1/3)*x^3+x-y^2);
                               /  1  3        2\
                            exp|- - x  + x - y |
                               \  3            /
> with(plots);
> plot3d(f, x = -2 .. 2, y = -1 .. 1, axes = box, scaling = constrained);

> dfdx := diff(f, x);
                       /  2    \    /  1  3        2\
                       \-x  + 1/ exp|- - x  + x - y |
                                    \  3            /
> dfdy := diff(f, y);
                                  /  1  3        2\
                          -2 y exp|- - x  + x - y |
                                  \  3            /
> evalf(solve({dfdx = 0, dfdy = 0}, {x, y}));
                     {y = 0., x = 1.}, {y = 0., x = -1.}
> fxx := diff(dfdx, x);
                                              2                     
                 /  1  3        2\   /  2    \     /  1  3        2\
         -2 x exp|- - x  + x - y | + \-x  + 1/  exp|- - x  + x - y |
                 \  3            /                 \  3            /
> fxy := diff(dfdx, y);
                        /  2    \      /  1  3        2\
                     -2 \-x  + 1/ y exp|- - x  + x - y |
                                       \  3            /
> fyx := diff(dfdy, x);
                        /  2    \      /  1  3        2\
                     -2 \-x  + 1/ y exp|- - x  + x - y |
                                       \  3            /
> fyy := diff(dfdy, y);
                   /  1  3        2\      2    /  1  3        2\
             -2 exp|- - x  + x - y | + 4 y  exp|- - x  + x - y |
                   \  3            /           \  3            /
> x := 1; y := 0;
                                      1
                                      0
> d := matrix([[fxx, fxy], [fyx, fyy]]);

> with(linalg);
> evalm(det(d));
                                           2
                                   /   /2\\ 
                                 4 |exp|-|| 
                                   \   \3// 
> x := -1; y = 0;
                                     -1
                                    0 = 0
> d := matrix([[fxx, fxy], [fyx, fyy]]);

> evalm(det(d));
                                            2
                                   /   /-2\\ 
                                -4 |exp|--|| 
                                   \   \3 // 
>