$
\def\l{\lambda} \def\o{{\tt1}}
\def\LR#1{\left(#1\right)}
\def\op#1{\operatorname{#1}}
\def\diag#1{\op{diag}\LR{#1}}
\def\Diag#1{\op{Diag}\LR{#1}}
\def\trace#1{\op{Tr}\LR{#1}}
\def\frob#1{\left\| #1 \right\|_F}
\def\qiq{\quad\implies\quad}
\def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}}
\def\c#1{\color{red}{#1}}
\def\CLR#1{\c{\LR{#1}}}
\def\fracLR#1#2{\LR{\frac{#1}{#2}}}
\def\gradLR#1#2{\LR{\grad{#1}{#2}}}
$First, let's give each variable a distinct name and remove those distracting subscripts (and denote row vectors by $v^T$ rather than $v$)
$$\eqalign{
x &= x_t^T,\quad &h=h_t^T,\quad &y = h_{t-1}^T,\quad &p = O_t^T \\
q &= b^T,\quad &R = W_h^T,\quad &S = W_x^T \\
}$$
Second, given a scalar function and its derivative
$\LR{f(z),\:f'=\grad fz}$ the differential is defined as
$$df = f'\,dz$$
These functions can be applied elementwise to a vector argument (like $p$) resulting in
$$h = f(p),\quad h'=f'(p),\qquad\qquad\quad dh = h'\odot dp$$
where $\odot$ denotes the elementwise/Hadamard product.
Such products can be replaced by diagonal matrices
$$H = \Diag{h},\quad H' = \Diag{h'},\qquad dh = H'\,dp$$
Fortunately, the derivative of tanh() is quite well known
$$dh = \LR{\o-h\odot h}\odot dp
\;=\; \CLR{I-H^2}dp \;\equiv\; \c{M}\,dp$$
We're almost done, we just need to calculate $dp$ and substitute
$$\eqalign{
p &= q + Ry + Sx \\
dp &= \LR{R\:dy+dR\:y} + \LR{S\:dx+dS\:x} \\
dh &= M\LR{R\:dy+dR\:y} + M\LR{S\:dx+dS\:x} \\
&= {MR\:dy+M\:dR\:y} + {MS\:dx+M\:dS\:x} \\
}$$
The vector gradients are easily identified from this expression
$$\eqalign{
\grad hy &= MR,\qquad \grad hx = MS \\
}$$
The gradients with respect to the matrix variables are trickier,
but they are
$$\eqalign{
\grad hR &= M\star y,\qquad \grad hS &= M\star x \\
}$$
where $\star$ denotes a dyadic/tensor product, i.e.
$$\eqalign{
\grad hS = M\star x \qiq \grad{h_i}{S_{jk}} = M_{ij}\,x_k \\
}$$
Such quantities are third-order tensors and are a PITA to work with using standard matrix-vector notation.
