This is a classical example of a Markov chain, a system evolving in time in discrete steps, each step depending only on the current value and not the previous values (ergodicity)
We will solve the problem as follows
Let $p_r^i$, $p_w^i$, $q_r^i$, $q_w^i$ be the expected quantities of red balls in first urn, white balls in first urn, red balls in second urn and white balls in second urn respectively. These may be fractional, as they are the product of quantity times probability. Let $i$ identify the time step. Given our initial conditions, at time step $i=0$
$p^0_r = 1$
$p^0_w = 2$
$q^0_r = 0$
$q^0_r = 5$
Now let us define a 4D vector, that glue together all these 4 quantities
$\vec{z}_i = \{p^i_r, q^i_r, p^i_w, q^i_w\}$
Note that I changed their order slightly for reasons that will become obvious later. What we will be trying to prove is that the vector $\vec{z}$ evolves linearly following the equation
$\vec{z}_{i+1} = M \vec{z}_{i}$
where $M$ is a fixed matrix. Once we have found this matrix, the problem can be solved by multiplying the initial vector by that matrix $n$ times
$\vec{z}_{n} = M^n \vec{z}_{0}$
and the final probability of the red ball being in the first urn being simply $p^n_r$
But let us now return to the definition of the matrix. Firstly, let us note that the total number of balls in each urn does not change
$n_p = p^i_r + p^i_w = 3$
$n_q = q^i_r + q^i_w = 5$
Let us consider the number of red balls in the first urn. It will decrease by 1 if a red ball is picked from this urn, and increased by 1 if a red ball is picked from another urn. We can write this as following
$p^{i+1}_r = p^{i}_r - \frac{1}{n_p} p^{i}_r + \frac{1}{n_q} q^i_r$
The third term is 1 times the probability of picking a red ball from the first urn, and the 4th term is 1 times the probability of picking a red ball from the second urn. By analogy
$p^{i+1}_r = p^{i}_r - \frac{1}{n_p} p^{i}_r + \frac{1}{n_q} q^i_r$
$p^{i+1}_w = p^{i}_w - \frac{1}{n_p} p^{i}_w + \frac{1}{n_q} q^i_w$
$q^{i+1}_r = q^{i}_r - \frac{1}{n_p} q^{i}_r + \frac{1}{n_q} p^i_r$
$q^{i+1}_w = q^{i}_w - \frac{1}{n_p} q^{i}_w + \frac{1}{n_q} p^i_w$
By reading off coefficients, we can now deduce the form of the matrix M
$M = \begin{bmatrix}
1 - 1/n_p & 1/n_q & 0 & 0 \\
1/n_p & 1 - 1/n_q & 0 & 0 \\
0 & 0 & 1 - 1/n_p & 1/n_q \\
0 & 0 & 1/n_p & 1 - 1/n_q \\
\end{bmatrix}$
As can be seen from the matrix, there is no coupling between red and white balls in terms of expected value evolution (off-diagonal blocks are zero), so the problem can be reduced to only considering red balls. Let us define a smaller 2D vector
$\vec{w}_i = \{p^i_r, q^i_r\}$
Then
$\vec{w}_{i+1} = N \vec{w}_{i}$
and
$N = \begin{bmatrix}
1 - 1/n_p & 1/n_q \\
1/n_p & 1 - 1/n_q \\
\end{bmatrix}$
Please tell me if you can figure it out from here. I can explain how one can accelerate the calculation of evolution by eigendecomposition of matrix N