How to solve the recurrence relation $a_{1}=2, a_{n}=\frac{a_{n-1}+2}{2 a_{n-1}+1}(n \geq 2)$ with generating functions? There's already a way to solve it, called "fixed point method", that is, from the relation we define its characteristic equation as $x=\dfrac{x+2}{2x+1}$，then we have $x_1=1,x_2=-1$. So the following relation established:
$$
\frac{a_{n}-1}{a_{n}+1}=\frac{\frac{a_{n-1}+2}{2 a_{n-1}+1}-1}{\frac{a_{n-1}+2}{2 a_{n-1}+1}+1}=-\frac{1}{3} \cdot \frac{a_{n-1}-1}{a_{n-1}+1}
$$
It is obvious that $\displaystyle \frac{a_{n}-1}{a_{n}+1}=\frac{1}{3} \cdot\left(-\frac{1}{3}\right)^{n-1}$, and then we have $a_{n}=\dfrac{3^{n}-(-1)^{n}}{3^{n}+(-1)^{n}}$.
My question is, how to solve this kind of recurrence relations with generating functions? Also, "fixed points" can be applied to solving recurrences like $a_{n+1}=\dfrac{a_{n}^{2}+b}{2 a_{n}+d}$, which seems impossible to solve using generating functions.

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\on}[1]{\operatorname{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
$\ds{\bbox[5px,#ffd]{%
\left\{\begin{array}{rcl}
\ds{a_{1}} & \ds{=} & \ds{2}
\\
\ds{a_{n}} & \ds{=} &
\ds{{a_{n - 1} + 2 \over 2 a_{n - 1} + 1}\,,\quad n \geq 2}
\end{array}\right.}}$

$\ds{\large Another\ Method}$: Lets $\ds{a_{n} \equiv x_{n}/y_{n}}$ such that
\begin{align}
&\bbox[5px,#ffd]{a_{n}  = {x_{n} \over y_{n}} =
{x_{n - 1}\,/y_{n - 1} + 2 \over 2x_{n - 1}\,/y_{n - 1} + 1} = {x_{n - 1} + 2y_{n - 1} \over 2x_{n - 1} + y_{n - 1}}}
\end{align}
and sets $\ds{\quad x_{n} = x_{n - 1} + 2y_{n - 1}\,,\quad
y_{n} = 2x_{n - 1} + y_{n - 1}}$ which is equivqlent to
\begin{align}
\pars{\begin{array}{c}\ds{x_{n}} \\ \ds{y_{n}}\end{array}}
& =
\pars{\begin{array}{cc}
\ds{1} & \ds{2}
\\
\ds{2} & \ds{1}
\end{array}}
\pars{\begin{array}{c}
\ds{x_{n - 1}} \\ \ds{y_{n - 1}}
\end{array}} =
\pars{\begin{array}{cc}
\ds{1} & \ds{2}
\\
\ds{2} & \ds{1}
\end{array}}^{2}
\pars{\begin{array}{c}
\ds{x_{n - 2}} \\ \ds{y_{n - 2}}
\end{array}}
\\[5mm] & = \cdots =
\pars{\begin{array}{cc}
\ds{1} & \ds{2}
\\
\ds{2} & \ds{1}
\end{array}}^{n - 1}
\pars{\begin{array}{c}
\ds{x_{1}} \\ \ds{y_{1}}
\end{array}} 
\end{align}
The above matrix has eigenvalues $\ds{\lambda_{1} = 3}$ and $\ds{\lambda_{2} = -1}$ with,
respectively, orthonormal eigenvectors
$\ds{{\bf u}_{1} = {1 \over \root{2}}{1 \choose 1}}$ and
$\ds{{\bf u}_{2} = {1 \over \root{2}}
{-1 \choose \phantom{-}1}}$. Then,
\begin{align}
\pars{\begin{array}{cc}
\ds{1} & \ds{2}
\\
\ds{2} & \ds{1}
\end{array}} & =
\sum_{j = 1}^{2}\lambda_{j}\,{\bf u}_{j}\,{\bf u}_{j}^{T}
\\[5mm] \mbox{and}\
\pars{\begin{array}{cc}
\ds{1} & \ds{2}
\\
\ds{2} & \ds{1}
\end{array}}^{n - 1} & =
\sum_{j = 1}^{2}\lambda_{j}^{n - 1}\,\,
{\bf u}_{j}\,{\bf u}_{j}^{T}
\\[2mm] & =
{ 1 \over 2}\pars{\begin{array}{cc}
\ds{3^{n -1} - \pars{-1}^{n}} & \ds{3^{n -1} + \pars{-1}^{n}}
\\
\ds{3^{n -1} + \pars{-1}^{n}} & \ds{3^{n -1} - \pars{-1}^{n}}
\end{array}}
\end{align}
Therefore,
\begin{align}
a_{n} & =
{\bracks{3^{n - 1} - \pars{-1}^{n}}\ \overbrace{x_{1}/y_{1}}^{\ds{= a_{1} = 2}}\ +\ 3^{n - 1} + \pars{-1}^{n} \over
\bracks{3^{n - 1} + \pars{-1}^{n}}x_{1}/y_{1} + 3^{n - 1} - \pars{-1}^{n}}
\\[5mm] & =
\bbx{3^{n} - \pars{-1}^{n} \over 3^{n} + \pars{-1}^{n}} \\ &
\end{align}
