I have a problem with the computation of the shape operator of a surface of revolution $f(u,v)=(\gamma_1(u) \cos(v), \gamma_2(u) \sin(v), \gamma_2(u))$, where $\gamma(u)=(\gamma_1(u),\gamma_2(u))$ is a an arclength-parametrized curve such that $\gamma_1(u)>0$.
I am using the following definitions
The shape operator is the smooth endomorphism field $A$ satisfying
$$ dN(X)=df(AX) $$
where $X \in TM$ with $TM$ denoting the tangent bundle of M. Here is my attempt at a solution. Let $p=(u,v) \in M$. The Jacobi matrix of $f$ at $p$ is given by
$$ f'(p)= \begin{pmatrix} \gamma_1'(u) \cos(v) & -\gamma_1(u) \sin(v) \\ \gamma_1'(u) \sin(v) & \gamma_1(u) \cos(v) \\ \gamma_2'(u) & 0 \end{pmatrix} $$
The unit normal is given by
$$ N=\frac{f_u \times f_v }{|f_u \times f_v |} =\frac{1}{\gamma_1(u)} \begin{pmatrix} -\gamma_2'(u) \gamma_1(u) \cos(v) \\ -\gamma_2'(u) \gamma_1(u) \sin(v) \\ \gamma_1'(u) \gamma_1(u) \end{pmatrix} =\begin{pmatrix} -\gamma_2'(u) \cos(v) \\ -\gamma_2'(u) \sin(v) \\ \gamma_1'(u) \end{pmatrix} $$
Then the Jacobi matrix of $N$ at $p$ is given by
$$ N'(p)= \begin{pmatrix} -\gamma_2''(u) \cos(v) & \gamma_2'(u) \sin(v) \\ -\gamma_2''(u) \sin(v) & -\gamma_2'(u) \cos(v) \\ \gamma_1''(u) & 0 \end{pmatrix} $$
Consider the coordinate vectorfields
$$ U: M \to TM , U(p)=(p,e_1) \\ V: M \to TM , V(p)=(p,e_2). $$.
Let $S=\begin{pmatrix} a & b \\ c & d \end{pmatrix}$ be the matrix representation of $A$. Then
$$ AU=aU+cV \\ AV=bU+dV. $$
Using the definition of $A$ we get $dN(U)=df(AU)$. With $dN(U)=N'(p) e_1, df(U)=f'(p)e_1, df(V)=f'(p) e_2$ and the linearity of $df$ we get
$$ \begin{pmatrix} -\gamma_2''(u) \cos(v) \\ -\gamma_2''(u) \sin(v) \\ \gamma_1''(u) \end{pmatrix} =a \begin{pmatrix} \gamma_1'(u) \cos(v) \\ \gamma_1'(u) \sin(v) \\ \gamma_2'(u) \end{pmatrix} + c \begin{pmatrix} -\gamma_1(u) \sin(v) \\ \gamma_1(u) \cos(v) \\ 0 \end{pmatrix} $$
Analogously, using $V$ instead of $U$ we get
$$ \begin{pmatrix} \gamma_2'(u) \sin(v) \\ -\gamma_2'(u) \cos(v) \\ 0 \end{pmatrix} =b \begin{pmatrix} \gamma_1'(u) \cos(v) \\ \gamma_1'(u) \sin(v) \\ \gamma_2'(u) \end{pmatrix} + d \begin{pmatrix} -\gamma_1(u) \sin(v) \\ \gamma_1(u) \cos(v) \\ 0 \end{pmatrix} $$
In total we get 6 equations
$$ -\gamma_2''(u) \cos(v)=a \gamma_1'(u) \cos(v)-c\gamma_1(u) \sin(v) \tag{1} $$
$$ -\gamma_2''(u) \sin(v)=a \gamma_1'(u) \sin(v)+c \gamma_1(u) \cos(v) \tag{2} $$
$$ \gamma_1''(u)=a \gamma_2'(u) \tag{3} $$
$$ \gamma_2'(u) \sin(v)=b \gamma_1'(u) \cos(v)-d \gamma_1(u) \sin(v) \tag{4} $$
$$ -\gamma_2'(u) \cos(v)=b \gamma_1'(u) \sin(v)+d\gamma_1(u) \cos(v) \tag{5} $$
$$ 0=b \gamma_2'(u) \tag{6} $$
Adding (1) to (2) yields
$$ -\gamma_2''(u) (\cos(v)+\sin(v))=a \gamma_1'(u) (\cos(v)+\sin(v))+c \gamma_1(u) (\cos(v)-\sin(v)). $$
This implies $c=0$. So we have the equation
$$ -\gamma_2''(u)=a \gamma_1'(u). \tag{7} $$
Dividing (1) by $\cos(v)$ and adding the resulting equation to (7) yields
$$ \gamma_1''(u)-\gamma_2''(u)=a(\gamma_1'(u)+\gamma_2'(u)). $$
Dividing by $(\gamma_1'(u)+\gamma_2'(u))$ and using the identity $1=\frac{(\gamma_1'(u)-\gamma_2'(u))}{(\gamma_1'(u)-\gamma_2'(u))}$ on the left-hand side yields
$$ a=\frac{-\gamma_2''(u) \gamma_1'(u)-\gamma_1''(u) \gamma_2'(u)}{(\gamma_1'(u))^2-(\gamma_2'(u))^2}. $$
Moreover (4) and (6) imply $b=0, d=-\frac{\gamma_2'(u)}{\gamma_1(u)}$. So the matrix representation of $A$ is given by
$$ S= \begin{pmatrix} \frac{-\gamma_2''(u) \gamma_1'(u)-\gamma_1''(u) \gamma_2'(u)}{(\gamma_1'(u))^2- (\gamma_2'(u))^2} & 0 \\ 0 & -\frac{\gamma_2'(u)}{\gamma_1(u)}. \end{pmatrix} $$
However according to this solution this result is not correct. According to this article I seem to have found the correct normal vector. This suggests that there must be an error in my subsequent computations but so far I have not been able to find any. Why does my approach not work?