I have the following two equations:
$$P_t = \frac{t-1}{t}P_{t-1} + \frac{1}{t}Q_{t-1}$$
$$Q_t = \frac{1}{t} + \frac{t-1}{t}Q_{t-1} - \frac{1}{t}P_{t-1}$$
with $P_0 = 0$ and $Q_0 = 0$. As time goes to infinity I know that both $P$ and $Q$ will equal 0.5. Notice that $t$ is time, so the system is not autonomous since time enters explicitly.
I have tried to solve it first by hand, exploiting the recursive structure, and even though there is a pattern at each step I cannot get an analytical expression for any $t$. I tried solving it on Matlab but have not succeeded: since the coefficients are time dependent, it seems I can only do this as a differential equation system and not as a difference equation system. Any suggestions? Is there a pen and paper method to solve this, or is there a way to do it Matlab or Python? Or can I cast this as a differential system?
This is a beautiful problem, as you are going to see in what follows.
First, you have correctly obtained the limit of the system, assuming that it exists. Now that we know it, let us rigorously show that it exists. Consider the change of unknowns $A_n = P_n - \frac 1 2, \ B_n = Q_n - \frac 1 2$ (motivated by the fact that $(\frac 1 2, \frac 1 2)$ is the limit of the system). We shall show that $A_n, B_n \to 0$.
Check that your system gets rewritten as
$$A_n = \left( 1 - \frac 1 n \right) A_{n-1} + \frac 1 n B_{n-1} \\ B_n = - \frac 1 n A_{n-1} + \left( 1 - \frac 1 n\right) B_{n-1} $$
with $A_0 = B_0 = - \frac 1 2$.
Rewrite this in matrix form as
$$\begin{pmatrix} A_n \\ B_n \end{pmatrix} = \begin{pmatrix} 1 - \frac 1 n & \frac 1 n \\ - \frac 1 n & 1 - \frac 1 n \end{pmatrix} \begin{pmatrix} A_{n-1} \\ B_{n-1} \end{pmatrix}$$
and use recursion to expand this as
$$\begin{pmatrix} A_n \\ B_n \end{pmatrix} = \begin{pmatrix} 1 - \frac 1 n & \frac 1 n \\ - \frac 1 n & 1 - \frac 1 n \end{pmatrix} \begin{pmatrix} 1 - \frac 1 {n-1} & \frac 1 {n-1} \\ - \frac 1 {n-1} & 1 - \frac 1 {n-1} \end{pmatrix} \dots \begin{pmatrix} 1 - \frac 1 1 & \frac 1 1 \\ - \frac 1 1 & 1 - \frac 1 1 \end{pmatrix} \begin{pmatrix} - \frac 1 2 \\ - \frac 1 2 \end{pmatrix} .$$
Note that we have been able to do this because the change of unknowns to $A_n$ and $B_n$ has annihilated the free term $\begin{pmatrix} 0 \\ \frac 1 n \end{pmatrix}$ that was present when working in matrix form with $P_n$ and $Q_n$.
Dont't even think to explicitly compute that product of matrices: first, because it is probably impossible; second, because you are going to see that it is not necessary, anyway.
Let $M_n = \begin{pmatrix} 1 - \frac 1 n & \frac 1 n \\ - \frac 1 n & 1 - \frac 1 n \end{pmatrix}$. Note that $M_n$ is of the form $\begin{pmatrix} a & b \\ - b & a \end{pmatrix}$ and that any two matrices of this form commute (check it!). Since all the matrices $(M_n)_{n \ge 1}$ commute, they may be simultaneously diagonalized over $\Bbb C$. This means that there exist an invertible matrix $X$ such that $M_n = X N_n X^{-1}$, where $N_n = \begin{pmatrix} x_n & 0 \\ 0 & \bar x_n \end{pmatrix}$ and $x_n, \bar x_n$ are the eigenvalues of $M_n$ (they are conjugate to each other because the characteristc equation of $M_n$ is of degree $2$ with real coefficients, therefore it has two complex-conjugated roots).
It is not difficult to compute the eigenvalues in the usual way and discover that $x_n = 1 - \frac 1 n + \frac {\rm i} n$. Why is this convenient? Because that ugly product of matrices becomes now a product of diagonal matrices, and this is easy to compute:
$$\begin{pmatrix} A_n \\ B_n \end{pmatrix} = X \begin{pmatrix} x_n & 0 \\ 0 & \bar x_n \end{pmatrix} \begin{pmatrix} x_{n-1} & 0 \\ 0 & \bar x_{n-1} \end{pmatrix} \dots \begin{pmatrix} x_1 & 0 \\ 0 & \bar x_1 \end{pmatrix} X^{-1} \begin{pmatrix} - \frac 1 2 \\ - \frac 1 2 \end{pmatrix} = \\ X \begin{pmatrix} x_n x_{n-1} \dots x_1 & 0 \\ 0 & \overline {x_n x_{n-1} \dots x_1} \end{pmatrix} X^{-1} \begin{pmatrix} - \frac 1 2 \\ - \frac 1 2 \end{pmatrix} .$$
Note that if we manage to show that $\prod \limits _{n=1} ^\infty x_n = \lim \limits _{n \to \infty} x_1 x_2 \dots x_n \to 0$ we are done. This is equivalent to showing that $\prod \limits _{n=1} ^\infty |x_n| = 0$. In order to study this, we shall take the logarithm of this product and show that $\sum \limits _{n=1} ^\infty \log |x_n| = - \infty$.
Using the value of $x_n$, the above series becomes $\sum \limits _{n=1} ^\infty \log \sqrt{\left( 1 - \frac 1 n \right)^2 + \frac 1 {n^2}}$; even simpler, since the square root gets out of the logarithm and of the series, let us study $\sum \limits _{n=1} ^\infty \log \left( 1 - \frac 2 n + \frac 2 {n^2} \right)$.
Note that this series has all the terms negative (because the arguments of $\log$ are numbers from $(0,1)$), so multiplying it by a minus sign turns it into a series with positive terms: $\sum \limits _{n=1} ^\infty - \log \left( 1 - \frac 2 n + \frac 2 {n^2} \right)$. The nice thing about this is that for series with positive terms we have the limit comparison convergence test: remembering that $\lim \limits _{x \to 0, x>0} \frac {- \log (1 - x)} x = 1$ and taking $x = \frac 2 n - \frac 2 {n^2}$, the limit comparison convergence test tells us that $\sum \limits _{n=1} ^\infty - \log \left( 1 - \frac 2 n + \frac 2 {n^2} \right)$ has the same behaviour as $\sum \limits _{n=1} ^\infty \left( \frac 2 n - \frac 2 {n^2} \right) = \sum \limits _{n=1} ^\infty \frac 2 n - \sum \limits _{n=1} ^\infty \frac 2 {n^2}$; the first term is clearly divergent summing up to $+ \infty$, and the second is convergent. This shows that $\sum \limits _{n=1} ^\infty \log \left( 1 - \frac 2 n + \frac 2 {n^2} \right) = - \infty$ (note the sign change, since we have reverted to our series with negative terms). Exponentiating it gives ${\rm e} ^{- \infty} = 0$, so $\prod \limits _{n=1} ^\infty x_n = 0$, which is what we were after.
To obtain explicit formulae for $A_n$ and $B_n$, let us go back to the diagonalization formula
$$\begin{pmatrix} A_n \\ B_n \end{pmatrix} = X \begin{pmatrix} x_n x_{n-1} \dots x_1 & 0 \\ 0 & \overline {x_n x_{n-1} \dots x_1} \end{pmatrix} X^{-1} \begin{pmatrix} - \frac 1 2 \\ - \frac 1 2 \end{pmatrix} .$$
Since $X$ and $X^{-1}$ diagonalize every matrix $M_n$, they diagonalize $M_1$, for which the computations are easy: the eigenvalues are ${\rm i}$ and $-{\rm i}$ and the corresponding eigenvectors are $\begin{pmatrix} 1 \\ {\rm i} \end{pmatrix}$ and $\begin{pmatrix} 1 \\ -{\rm i} \end{pmatrix}$, therefore $M_n = X \begin{pmatrix} x_n & 0 \\ 0 & \bar x_n \end{pmatrix} X^{-1}$, where $X = \begin{pmatrix} 1 & 1 \\ {\rm i} & -{\rm i} \end{pmatrix}$ and $X^{-1} = \frac 1 2 \begin{pmatrix} 1 & -{\rm i} \\ 1 & {\rm i} \end{pmatrix}$. Plugging these ingredients into the above expression, we get:
$$\begin{pmatrix} A_n \\ B_n \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ {\rm i} & -{\rm i} \end{pmatrix} \begin{pmatrix} x_n x_{n-1} \dots x_1 & 0 \\ 0 & \overline {x_n x_{n-1} \dots x_1} \end{pmatrix} \frac 1 2 \begin{pmatrix} 1 & -{\rm i} \\ 1 & {\rm i} \end{pmatrix} \begin{pmatrix} - \frac 1 2 \\ - \frac 1 2 \end{pmatrix} = - \frac 1 2 \begin{pmatrix} \text{Re } x_1 x_2 \dots x_n + \text{Im } x_1 x_2 \dots x_n \\ \text{Re } x_1 x_2 \dots x_n - \text{Im } x_1 x_2 \dots x_n \end{pmatrix} ,$$
whence it finally follows that
$$P_n = \frac {1 - \text{Re } \prod \limits _{k=1} ^n \left( 1 - \frac 1 k + \frac {\rm i} k \right) - \text{Im } \prod \limits _{k=1} ^n \left( 1 - \frac 1 k + \frac {\rm i} k \right)} 2, \quad Q_n = \frac {1 - \text{Re } \prod \limits _{k=1} ^n \left( 1 - \frac 1 k + \frac {\rm i} k \right) + \text{Im } \prod \limits _{k=1} ^n \left( 1 - \frac 1 k + \frac {\rm i} k \right)} 2 .$$
I strongly doubt that one can find a closed formula for these.