There probably is no general algorithm, but one could say that the following approach is rather general:

• Rewrite the function $F\colon \mathbb{R}^{2}\to\mathbb{R}^{2}$ as an affine map: $$ F\left(\begin{pmatrix}a\\ b\end{pmatrix} \right) = \begin{pmatrix}1 - br + a\\ 1 - br + a + b\end{pmatrix} = \underbrace{\begin{pmatrix}1\\ 1\end{pmatrix}}_{=: f} + \underbrace{\begin{pmatrix}1 & -r\\ 1 & 1-r\end{pmatrix}}_{=: A}. \begin{pmatrix}a\\ b\end{pmatrix} $$
• You the ansatz that $G\colon \mathbb{R}^{2}\to\mathbb{R}^{2}$ is also affine, i.e. $G(v) = g + Bv$ for some $g\in\mathbb{R}^{2}$ and $B\in \mathbb{R}^{2\times 2}$ and compute $$ G(G(v)) = g + B G(v) = g + Bg + B^2v $$
• Comparing this with the above formula you arrive at the conditions $$ g + Bg \stackrel{!}{=} f, \qquad B^2 = A. $$
• Solve the second equation first to obtain $B$ and then plug the result into the first equation to obtain $g = (\mathrm{Id}+B)^{-1}f$ and you are done.

However, solving $B^2 = A$ might not be easy in general. In your case you could diagonalize $A$ first, $$ A = S^{-1} \begin{pmatrix}\lambda_{1} & 0 \\ 0 & \lambda_2\end{pmatrix} S, $$ but note that, depending on the sign of $r$, the eigenvalues $\lambda_{j}$ of $A$ might be complex. Then choose $$ B = S^{-1} \begin{pmatrix}\sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_2}\end{pmatrix} S $$ which obviously solves $B^2 = A$. It does not matter which square roots you choose, but $B$ might be a complex matrix. I leave the details to you.

There probably is no general algorithm, but one could say that the following approach is rather general:

• Rewrite the function $F\colon \mathbb{R}^{2}\to\mathbb{R}^{2}$ as an affine map: $$ F\left(\begin{pmatrix}a\\ b\end{pmatrix} \right) = \begin{pmatrix}1 - br + a\\ 1 - br + a + b\end{pmatrix} = \underbrace{\begin{pmatrix}1\\ 1\end{pmatrix}}_{=: f} + \underbrace{\begin{pmatrix}1 & -r\\ 1 & 1-r\end{pmatrix}}_{=: A}. \begin{pmatrix}a\\ b\end{pmatrix} $$
• You the ansatz that $G\colon \mathbb{R}^{2}\to\mathbb{R}^{2}$ is also affine, i.e. $G(v) = g + Bv$ for some vector $g$ and matrix $B$ and compute $$ G(G(v)) = g + B G(v) = g + Bg + B^2v $$
• Comparing this with the above formula, your requirement that $F(v) = G(G(v))$ yields the conditions $$ g + Bg \stackrel{!}{=} f, \qquad B^2 = A. $$
• Solve the second equation first to obtain $B$ and then plug the result into the first equation to obtain $g = (\mathrm{Id}+B)^{-1}f$ and you are done.

However, solving $B^2 = A$ might not be easy in general. In your case you could diagonalize $A$ first, $$ A = S^{-1} \begin{pmatrix}\lambda_{1} & 0 \\ 0 & \lambda_2\end{pmatrix} S, \qquad S \in \mathbb{C}^{2\times 2}, $$ but note that, depending on the sign of $r$, the eigenvalues $\lambda_{j}$ of $A$ might be complex. Then choose $$ B = S^{-1} \begin{pmatrix}\sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_2}\end{pmatrix} S $$ which obviously solves $B^2 = A$. It does not matter which square roots you choose, but $B$ might be a complex matrix. I leave the details to you.

Remark: For $r=0$, $A$ is not diagonalizable, but it is easy to see that in this case $$ A = \begin{pmatrix}1 & 0\\ 1 & 1\end{pmatrix} = B^2 \quad\text{for}\quad B = \begin{pmatrix}1 & 0\\ 1/2 & 1\end{pmatrix}.$$