What's wrong with my understanding of how RNNs work?

Recently, I've been trying to derive the mathematics behind various Neural Network structures. I managed to derive the MLP and tested it to be on par with a Keras implementation (Using the MNIST dataset).

Moving on to RNNs, I assume that the only additional mathematics required is for the recurrent weight matrix $$R$$.

I take the convention $$x_{i}^{t}$$ to be the activation of the $$i$$th layer of neurons at time step $$t$$.

The following mathematics describes my network. Specifically, it describes the feedforward process.

$$x_{i}^{t}=f(x_{i-1}^{t}W_{i-1}+b_{i-1}+x_{i}^{t-1}R_{i})$$ $$z_{i}^{t}=x_{i-1}^{t}W_{i-1}+b_{i-1}+x_{i}^{t-1}R_{i},$$

where $$f$$ is some activation function.

For backpropagation, the gradients of $$W_{i}$$ and $$b_{i}$$ are the same as for MLP. There is only a need to calculate $$\frac{\partial L}{\partial R_{i}}$$, where $$L$$ is some loss function. For clarity, I did everything index-wise.

$$\frac{\partial L}{\partial R_{i}}=\begin{pmatrix} \frac{\partial L}{\partial R_{i}[0,0]} & \frac{\partial L}{\partial R_{i}[0,1]} & \dots\\ \frac{\partial L}{\partial R_{i}[1,0]} & \ddots\\ \vdots \end{pmatrix}$$

$$\frac{\partial L}{\partial R_{i}[j,k]}=\sum_{n}\left(\frac{\partial L}{\partial x_{i}^{t}[n]}\frac{\partial x_{i}^{t}[n]}{\partial R_{i}[j,k]}\right)$$

$$\frac{\partial L}{\partial x_{i}^{t}}$$ should have been calculated via the standard MLP procedure. That leaves $$\frac{\partial x_{i}^{t}[n]}{\partial R_{i}[j,k]}$$

$$x_{i}^{t}[n]=f(x_{i-1}^{t}W_{i-1}[:,n]+b_{i-1}[n]+x_{i}^{t-1}R_{i}[:,n])$$

$$\frac{\partial x_{i}^{t}[n]}{\partial R_{i}[j,k]}=f'(z_{i}^{t}[n])\left[\frac{\partial x_{i}^{t-1}}{\partial R_{i}[j,k]}R_{i}[:,n]+\delta_{n,k}x_{i}^{t-1}[j]\right]$$ Here, $$\frac{\partial x_{i}^{t-1}}{\partial R_{i}[j,k]}$$ should have been calculated at the time step before. It is the source of recurrence and is the term that blows up or decays to cause exploding/vanishing gradients. The training set is then fed into the network one by one and gradients are updated at every time step. I was fairly confident in this derivation but no one online seems to derive BPTT this way. I tried implementing this and testing it on the XOR training set like the original paper by Elman did, but it showed no improvement over just using my MLP. Additionally, a Keras implementation of a similar architecture handily beats my own implementation.

Is there something wrong with my derivation? It seems pretty intuitive to me but I find it very unsettling that I can't find another person online who does it this way. Or perhaps there's something fundamentally wrong about my understanding of how RNNs work?