Suppose $M \subset \mathbb R^{n+1}$ is the graph of a function $f : U \to \mathbb R$, for some open subset $U \subset \mathbb R^n$. I'm trying to find the coefficients of the matrix of the shape operator $s : T_p M \to T_p M$ (for $p \in M$) in terms of graph coordinates. These are global coordinates of $M$ given by $\phi : M \to U$, $\phi(u, f(u)) = u$.
What I've tried. Letting $(x^1, \ldots, x^{n+1})$ denote the usual Cartesian coordinates in $\mathbb R^{n+1}$, the coordinate tangent vectors of $M$ in graph coordinates are related to the coordinate tangent vectors of $\mathbb R^{n+1}$ via $$ \frac{\partial}{\partial u^j} = \frac{\partial}{\partial x^j} + \frac{\partial f}{\partial x^j} \frac{\partial}{\partial x^{n+1}} $$ (This can be seen by considering the map $\phi^{-1} : U \to M$, $\phi^{-1}(u) = (u, f(u))$ and explicitly calculating $\frac{\partial\left(\phi^{-1}\right)^j}{\partial u^i} \frac{\partial}{\partial x^j}$.) The upward unit normal field of $M$ is given by $$ N = \frac{\mathrm{grad} F}{|\mathrm{grad} F|}, $$ where $F : \mathbb R^{n+1} \to \mathbb R$ is the function $F(x,x^{n+1}) = x^{n+1} - f(x)$ (so $M$ is the level set of $F$ for $0$). Therefore, $$ N = \left(1+\sum_{i=1}^n \left(\frac{\partial f}{\partial x^i}\right)^2\right)^{-1/2}\left(-\sum_{i=1}^n \frac{\partial f}{\partial x^i}\frac{\partial}{\partial x^i} + \frac{\partial}{\partial x^{n+1}}\right) = -\sum_{i=1}^n \frac{\partial f/\partial x^i}{|\mathrm{grad} F|} \frac{\partial}{\partial x^i} + \frac 1{|\mathrm{grad} F|} \frac{\partial}{\partial x^{n+1}} $$ In particular, if we denote $N = N^i \frac{\partial}{\partial x^i}$, then $\frac{\partial}{\partial x^{n+1}} N^i \equiv 0$ for every $i$, and thus $\overline\nabla_{\partial/\partial x^{n+1}} N = 0$ (where $\overline\nabla$ is the Levi-Civita connection of the Euclidean metric on $\mathbb R^{n+1}$). Meanwhile, by direct computation, for $1 \leq i, j \leq n$, $$ \frac{\partial}{\partial x^j} N^i = -\frac 1{|\mathrm{grad} F|^3} \frac{\partial f}{\partial x^i}\sum_{k=1}^n \frac{\partial f}{\partial x^k} \frac{\partial^2f}{\partial x^k \partial x^j} + \frac{1}{|\mathrm{grad} F|} \frac{\partial^2 f}{\partial x^i \partial x^j} $$ and $$ \frac{\partial}{\partial x^j} N^{n+1} = -\frac 1{|\mathrm{grad} F|^3}\sum_{i=1}^n\frac{\partial f}{\partial x^i}\frac{\partial^2 f}{\partial x^i \partial x^j} $$ The Weingarten equation for hypersurfaces says $sX = -\overline\nabla_X N$ for $X \in \mathcal X(M)$. So using what we just calculated, \begin{align} s\frac{\partial}{\partial u^j} &= -\overline\nabla_{\partial/\partial x^j} N - \frac{\partial f}{\partial x^j} \overline\nabla_{\partial/\partial x^{n+1}} N = -\overline\nabla_{\partial/\partial x^j} N = -\sum_{i=1}^n \left(\frac{\partial}{\partial x^j} N^i\right) \frac{\partial}{\partial x^i}\\ &= \sum_{i=1}^n\left(\frac 1{|\mathrm{grad} F|^3} \frac{\partial f}{\partial x^i}\sum_{k=1}^n \frac{\partial f}{\partial x^k} \frac{\partial^2f}{\partial x^k \partial x^j} - \frac{1}{|\mathrm{grad} F|} \frac{\partial^2 f}{\partial x^i \partial x^j}\right)\frac{\partial}{\partial x^i} \\ &\qquad \qquad \qquad + \frac 1{|\mathrm{grad} F|^3} \sum_{i=1}^n \frac{\partial^2 f}{\partial x^i \partial x^j} \frac{\partial f}{\partial x^i} \frac{\partial}{\partial x^{n+1}} \\ &= \frac 1{|\mathrm{grad} F|^3} \sum_{i=1}^n \left[ \left(\frac{\partial f}{\partial x^i}\sum_{k=1}^n \frac{\partial f}{\partial x^k} \frac{\partial^2f}{\partial x^k \partial x^j} - |\mathrm{grad}F|^2 \frac{\partial^2 f}{\partial x^i \partial x^j}\right)\frac{\partial}{\partial x^i} +\frac{\partial^2 f}{\partial x^i \partial x^j} \frac{\partial f}{\partial x^i} \frac{\partial}{\partial x^{n+1}} \right] \end{align} My problem: In principle, I should be able to now express $s\frac{\partial}{\partial u^i}$ in terms of the basis of coordinate vector fields given by $\frac{\partial}{\partial u^j} = \frac{\partial}{\partial x^j} + \frac{\partial f}{\partial x^j} \frac{\partial}{\partial x^{n+1}}$. This will give me the coefficients for $s_{ij}$ I need. But it's not clear to me that I'm able to do that given what I just found. On the other hand, I'm reasonably sure that the expression I've found is orthogonal to $N$, so it lies in the tangent space of $M$ at each point $p \in M$. Am I going wrong somewhere? Is there a better way to find the coefficients of the shape operator?
EDIT: I now see I was trying to plug a round peg into a square hole. My computations were correct, and agree with the answer below (though are written less elegantly). The coefficients in graph coordinates are $$ s_{ij} = \frac 1{|\mathrm{grad} F|^3} \frac{\partial f}{\partial x^i}\sum_{k=1}^n \frac{\partial f}{\partial x^k} \frac{\partial^2f}{\partial x^k \partial x^j} - \frac{1}{|\mathrm{grad} F|} \frac{\partial^2 f}{\partial x^i \partial x^j} $$ since $\frac{\partial}{\partial u^j} = \frac{\partial}{\partial x^j} + \frac{\partial f}{\partial x^j} \frac{\partial}{\partial x^{n+1}}$. In particular, comparing this coordinate expression to the $\frac{\partial}{\partial x^{n+1}}$ term above implies $$\sum_{i=1}^n s_{ij} \frac{\partial f}{\partial x^i} = \frac 1{|\mathrm{grad}F|^3} \sum_{k=1}^n \frac{\partial f}{\partial k} \frac{\partial^2 f}{\partial x^k \partial x^j}$$ although explicitly writing out why the above is true would be rather challenging.
If $X(u) = (u,f(u))$ is the usual Monge parametrization for $M = {\rm gr}(f)$, we have that $$DX(u) = ({\rm Id}, Df(u)),$$and thus a unit vector field in $\Bbb R^{n+1}$ normal to $M$ is just $$N(u) = \frac{(-\nabla f(u), 1)}{\sqrt{1+\|\nabla f(u)\|^2}}.$$The shape operator is given by $s = -{\rm d}N$. But to find its components, we may start with the scalar version of the second fundamental form, whose coefficients are $h_{ij} = \langle X_{ij}, N\rangle$. But if $(e_1,\ldots, e_n)$ is the standard basis for $\Bbb R^n$, we have that $$X_{ij}(u) = D^2X(u)(e_i,e_j)= (0, D^2f(u)(e_i,e_j)),$$and so $$h_{ij}(u) =\left\langle(0, D^2f(u)(e_i,e_j)), \frac{(-\nabla f(u), 1)}{\sqrt{1+\|\nabla f(u)\|^2}} \right\rangle = \frac{D^2f(u)(e_i,e_j)}{\sqrt{1+\|\nabla f(u)\|^2}}.$$Thus $${\rm II}_{X(u)} = \frac{D^2f(u)}{\sqrt{1+\|\nabla f(u)\|^2}},$$where $D^2f(u)$ is the Hessian bilinear form of $f$ at $u$. Now, we have that $$s(X_j) = \sum_i h^i_{~j}X_i,$$where $h^i_{~j} = \sum_k g^{ik}h_{kj}$. But $$g_{ij} = \delta_{ij} + (\partial_if)(\partial_jf) \implies g^{ij} = \delta^{ij} - \frac{(\partial_if)(\partial_jf)}{1+\|\nabla f\|^2}.$$With this we compute $$h^i_{~j} = \sum_k \left(\delta^{ik} - \frac{(\partial_if)(\partial_kf)}{1+\|\nabla f\|^2}\right) \frac{\partial_k\partial_jf}{\sqrt{1+\|\nabla f\|^2}} = \frac{\partial_i\partial_jf}{\sqrt{1+\|\nabla f\|^2}} - \sum_k \frac{(\partial_if)(\partial_kf)(\partial_k\partial_jf)}{(1+\|\nabla f\|^2)^{3/2}}.$$A bit more of work shows that $H = 0$ if and only if $${\rm div}\left(\frac{\nabla f}{\sqrt{1+\|\nabla f\|^2}}\right) = 0,$$which is the minimal hypersurface equation in divergence form.