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}.$$