Introducir secuencias de $A_n, B_n, C_n$, cuya escala general a ser determinado, de tal manera que
$a_n = A_n/C_n$, $b_n = B_n/C_n$. Podemos reescribir nuestra ecuaciones como
$$
\requieren{cancel}
\begin{cases}
a_{n+1} &= \frac{a_n+1}{a_n+b_n+3}\\
b_{n+1} &= \frac{b_n+1}{a_n+b_n+3}
\end{casos}
\quad\ffi\quad
\begin{pmatrix}A\\B\\C\end{pmatrix}_{n+1}
= \color{red}{\cancelto{1}{\color{gris claro}{\frac{C_{n+1}}{A_n+B_n+3C_n}}}}
\begin{pmatrix}1&0&1 \\ 0&1&1 \\ 1&1&3 \end{pmatrix}
\begin{pmatrix}A\\B\\C\end{pmatrix}_{n}
$$
Si elegimos la escala para hacer $C_{n+1} = A_n+B_n+3C_n$ e imponer la condición inicial
$(A_1,B_1,C_1) = (2,1,1)$ sobre ellos, el correspondiente $a_n = A_n/C_n$ $b_n = B_n/C_n$ va a resolver el problema original.
Deje $M$ $3\times 3$ matriz aparecen arriba. Es fácil de comprobar
$$
M\begin{pmatrix}1 \\ -1 \\ 0\end{pmatrix} = \begin{pmatrix}1 \\ -1 \\ 0\end{pmatrix}
\quad\text{ y }\quad
M\begin{pmatrix}1 \\ 1 \\ \rho\end{pmatrix} = (1+\rho) \begin{pmatrix}1 \\ 1 \\ \rho\end{pmatrix}
\quad\text{ para }\rho = 1\pm \sqrt{3}.
$$
Vamos $\alpha = 1+\sqrt{3}$, $\beta = 1-\sqrt{3}$, podemos expresar el vector inicial como
$$
\begin{pmatrix} A \\ B \\ C\end{pmatrix}_1 =
\begin{pmatrix} 2 \\ 1 \\ 1\end{pmatrix} =
\frac12\begin{pmatrix}1 \\ -1 \\ 0\end{pmatrix}
+ \frac{3+\alpha}{2\sqrt{3}\alfa} \begin{pmatrix}1 \\ 1 \\ \alpha\end{pmatrix}
- \frac{3+\beta}{2\sqrt{3}\beta} \begin{pmatrix}1 \\ 1 \\ \beta\end{pmatrix}
$$
Aviso $\alpha^2 = 2(1+\alpha)$, $\beta^2 = 2(1+\beta)$, nos encontramos
$$
\begin{pmatrix} A \\ B \\ C\end{pmatrix}_{n+1} =
\frac12\begin{pmatrix}1 \\ -1 \\ 0\end{pmatrix}
+ \frac{1}{2^{n+1}\sqrt{3}}\left\{
(3+\alpha)\alpha^{2n-1} \begin{pmatrix}1 \\ 1 \\ \alpha\end{pmatrix}
- (3+\beta)\beta^{2n-1} \begin{pmatrix}1 \\ 1 \\ \beta\end{pmatrix}
\right\}
$$
De estos, encontramos los términos queremos suma puede ser transformado como
$$\begin{align}
\frac{(a_n+1)(b_n+1)}{a_n+b_n+3}
= & \frac{A_{n+1}B_{n+1}}{C_{n+1}C_n}\\
= & \frac{\left((3+\alpha)\alpha^{2n-1} - (3+\beta)\beta^{2n-1}\right)^2
- 3\cdot 4^n}{
2\left((3+\alpha)\alpha^{2n} - (3+\beta)\beta^{2n}\right)
\left((3+\alpha)\alpha^{2n-2} - (3+\beta)\beta^{2n-2}\right)
}
\end{align}
$$
Vamos $\lambda = \frac{\beta}{\alpha}$, $\mu = \frac{3+\beta}{3+\alpha}$ y aviso a $\alpha\beta = -2$, se puede simplificar por encima de desastre como
$$\begin{align}
\frac{\left(1 -\mu\lambda^{2n-1}\right)^2 + \frac{6}{(3+\alpha)^2}\lambda^{2n-1}}
{2\left(1 -\mu\lambda^{2n}\right)\left(1 -\mu\lambda^{2n-2}\right)}
= & \frac12 + \frac{\Delta \lambda^{2n-1}}{(1 -\mu\lambda^{2n})(1 -\mu\lambda^{2n-2})}\\
= & \frac12 + \Delta' \left\{\frac{1}{1 -\mu\lambda^{2n}} - \frac{1}{1 -\mu\lambda^{2n-2}}\right\}
\end{align}$$
donde
$$2\Delta = \mu(\lambda+\lambda^{-1}-2) + \frac{6}{(3+\alpha)^2}
\quad\text{ y }\quad
\Delta' = \frac{\Delta}{\mu(\lambda\lambda^{-1})} = -\frac{6\sqrt{3}}{13}.$$
Por tanto, tenemos
$$\begin{align}
\sum_{k=1}^n\frac{(a_k+1)(b_k+1)}{a_k+b_k+3}
= & \frac{n}{2} + \frac{6\sqrt{3}}{13}\left(\frac{1}{1-\mu}-\frac{1}{1-\mu\lambda^{2n}}\right)\\
= & \frac{n}{2} + \frac{6\sqrt{3}}{13}\left(\frac{1}{1-\frac{4-\sqrt{3}}{4+\sqrt{3}}}-\frac{1}{1-\frac{4-\sqrt{3}}{4+\sqrt{3}}(2-\sqrt{3})^{2n}}\right)\\
= & \frac{n}{2} + \frac{3}{4+\sqrt{3}}\left( 1 - \frac{2\sqrt{3}}{(4+\sqrt{3})(2+\sqrt{3})^{2n} - (4-\sqrt{3})}\right)
\end{align}$$