Limit of probability distribution of $X_{i+1},Y_{i+1} \overset{d}= f(X_i,Y_i)$ for $i \rightarrow \infty$ The Question
Imagine a setup where you start with two iid random variables $X_0,Y_0\sim\mathcal D_0$, where $\mathcal D_0$ is some discrete distribution over $[N]$ for some $N\in\mathbb N$. You also have a function $f:N\times N\longrightarrow N$, and you define $\mathcal D_{i+1}$ to be the distribution of $f(X_i, Y_i)$, and take $X_{i+1},Y_{i+1} \sim \mathcal D_{i+1}$. For which $f$ does the limit $\lim_{i\rightarrow\infty} \mathcal D_i$ exist? How to calculate it?
For the sake of simplicity you can assume that $\mathcal D_0$ is the uniform distribution over $[N]$.
An Attempt
I thought of "pivoting" the problem, in the sense of defining $p_i(n) := \mathrm{Prob}( X_i = n)$, which means
$$
\begin{split}
    p_{i+1}(n) = \mathrm{Prob}(f(X_i,X'_i)=n) &= \!\!\!\!\!\!\sum_{(m,k) \in f^{-1}(n)}\!\!\!\!\!\! p_i(m)p_i(k) \\
 &= \sum_{m\in[N]}\sum_{k\in[N]} p_i(m) p_i(k)
                 \begin{cases}
                      1 & n=f(m,k) \\ 0 & \text{otherwise.}
                 \end{cases}
\end{split}
$$
Maybe this set of coupled recurrence relations can somehow be translated into a set of ODEs, and then one could apply Picard-Lindelöf?
I'd be happy for any pointers!!!
PS: I heavily edited the question after realising that there was no good way of formulating it with balls and colors. Thanks for pointing that out!
 A: Not an answer / too long for a comment.
Your equations are necessary for any limits that exist, but they are not sufficient to show limits exist.  This is the difference between a stationary solution vs actually tending to such a solution from a given starting point.
Here is a simple $f$ which does not even depend on its second argument, where the limit does not exist.  Consider $N \ge 3$.  For any $y$:

*

*$f(1, y) = N$


*$f(k, y) = 1$ where $k \neq 1$.
Starting with the uniform distribution, all future distributions will have only $p(1)$ and $p(N)$ be non-zero, but their values will oscillate between $\frac1N$ and $\frac{N-1}{N}$.  The stationary solution is of course both values equal $\frac12$ but this will not be reached from the uniform starting distribution.
Obviously my answer is inspired by Markov Chains, where stationary and limit are related but not the same thing.  So if you are looking for limiting behavior, there are other things to consider besides the "quadratic map" issue.
