I'll provide a sketch of the derivation. Omitting the bias term (since anyways we take derivatives later), the recursion looks like:
$$\mathbf{h_{t+1}}=tanh(\mathbf{ U x_{t}+W h_{t}}) $$ where the $tanh$ is taken elementwise.

Now, since $\mathbf{h_{t}}$ and $\mathbf{h_{t+1}}$ are vectors, the derivative $\frac{\partial \mathbf{h_{t+1}}}{\partial \mathbf{h_{t}}}$ is a  [Jacobian][1]. 

 1. Notice that if $y=tanh(x), dy/dx=1-tanh^2(x)=1-y^2$.
 2. Let's see how a single element of the Jacobian looks like. Assume that the hidden layers are of dimension $n$ Now $h_{t+1, i}=tanh(\sum_{k=1} ^{n} w_{ik}h_{t,k}+ g(x))$. Here $g(x)$ stands for some function of x.
 
 3. Hence $\frac{\partial h_{t+1, i}}{\partial h_{t, j}}$ = $1-h_{t+1,i}^{2}w_{i,j}$. This corresponds to element in position $(i,j)$ of the Jacobian (from #1 above).
 4. This multiplier $1-h_{t+1,i}^{2}$ applies to each element of the $i^{th}$ row of the Jacobian.
 5.  From basic linear algebra you can show that pre-multiplying a matrix $A$ by a diagonal matrix $diag(b_{i,i})$ is equivalent to scaling up row $i$ of matrix $A$ by $b_{i,i}$

**Hence the term $diag(1-\mathbf{h_{t+1}^{2}})$ occurs as a premultiplier to the matrix $\mathbf{W}$**


  [1]: https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant