This seems trivial but I am stuck. To gain intuition for a bigger problem, I am trying to compute the shape operator of the graph of $f(x,y)=xy$ at the point $p=(0,0,0)$, call this surface $\Sigma$. In class we defined the shape operator at $p\in \Sigma$, in terms of the covariant derivative.
$$S(\vec v) = -\nabla_{\vec v}N=-\frac{d}{dt}\bigg|_{t=0}N(p+t\vec v).$$ where $\vec v\in T_p\Sigma$ and $N$ is normal to $\Sigma$.
Let $F(x,y,z)=xy-z$. I know $N = \frac{\nabla F}{|\nabla F|}=\frac{1}{\sqrt{1+x^2+y^2}}(y,x,-1)$. Say $\vec e_1,\vec e_2$ are basis vectors for $\Sigma$. I know I need to find $S(\vec e_1)$ and $S(\vec e_2)$, from there I can find the matrix of $S$. This is where I am stuck, if someone can show how to do the calculation, that would be helpful.
The surface is naturally parameterized by $(x, y)$: $$(x, y)\mapsto (x, y, xy)$$ So there are natural tangent vectors $$e_1=\frac{\partial}{\partial x}=(1, 0, y)$$ $$e_2=\frac{\partial}{\partial y}=(0, 1, x)$$ Hence (using your notation and let $r=\sqrt{1+x^2+y^2}$) $$\begin{align*} S(e_1) &=-\frac{d}{dt} N(p+te_1)\\ &= -\frac{\partial N}{\partial x}\\ &= -(-xy, 1+y^2, x)/r^3\\ S(e_2) &=-\frac{d}{dt} N(p+te_2)\\ &= -\frac{\partial N}{\partial y}\\ &= -(1+x^2, -xy, y)/r^3 \end{align*}$$ Evaluating at $x=y=0$ we get $$S(e_1)|_{(0, 0)}=-(0, 1, 0)=-e_2, S(e_2)_{(0, 0)}=-(1, 0, 0)=-e_1$$ Hence the matrix is $$\begin{bmatrix}0 &-1 \\ -1 & 0\end{bmatrix}$$