I'll provide a sketch of the derivation. Omitting the bias term (since anyways we take derivatives later), the recursion looks like: $$h_{t+1}=tanh(x_{t} U+W h_{t}) =tanh(p), \textrm{say}$$$$\mathbf{h_{t+1}}=tanh(\mathbf{ U x_{t}+W h_{t}}) $$ where the tanh$tanh$ is taken elementwise.
Now, since $h_{t}$$\mathbf{h_{t}}$ and $h_{t+1}$$\mathbf{h_{t+1}}$ are vectors, the derivative $\frac{\partial h_{t+1}}{\partial h_{t}}$$\frac{\partial \mathbf{h_{t+1}}}{\partial \mathbf{h_{t}}}$ is a Jacobian. Now,
Notice that if $y=tanh(x), dy/dx=1-tanh^2(x)=1-y^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.
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).
This multiplier $1-h_{t+1,i}^{2}$ applies to each element of the $i^{th}$ row of the Jacobian.
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}$
Hence the term $diag(1-\mathbf{h_{t+1}^{2}})$ occurs as a premultiplier to the matrix $\mathbf{W}$