I think you might be referring to what is also covered in the paper referenced by the documentation of the matlab function dare. In that paper they call it the Schur method. This is based on adding a co-state, similar to Pontryagin's maximum (or minimum) principle.
The solution to the discrete algebraic Riccati equation (DARE) is used in the infinite horizon discrete LQR problem, which can be formulated as follows
$$
\begin{align}
\min_u & \frac{1}{2} \sum_{k=0}^\infty x_k^\top Q\,x_k + u_k^\top R\,u_k \\
\text{s.t.}\, &\, x_{k+1} = A\,x_k + B\,u_k
\end{align} \tag{1}
$$
By introducing Lagrange multipliers this can also be written as
$$
\min_u J = \sum_{k=0}^\infty \frac{1}{2} \left(x_k^\top Q\,x_k + u_k^\top R\,u_k\right) + \lambda_k^\top(A\,x_k + B\,u_k - x_{k+1}) \tag{2}
$$
The optimal solution can be found by setting the partial derivative of the cost function $J$ with respect to each variable to zero
$$
\left.
\begin{align}
\frac{\partial\,J}{\partial\,x_k} &= Q\,x_k + A^\top \lambda_k - \lambda_{k-1} = 0 \\
\frac{\partial\,J}{\partial\,\lambda_k} &= A\,x_k + B\,u_k - x_{k+1} = 0 \\
\frac{\partial\,J}{\partial\,u_k} &= R\,u_k + B^\top \lambda_k = 0
\end{align}
\right\} \forall\, k \geq 0 \tag{3}
$$
Solving the last for $u_k$ gives
$$
u_k = -R^{-1} B^\top \lambda_k. \tag{4}
$$
By substituting $(4)$ into the remaining equations of $(3)$ it can be shown those can be written as
$$
\begin{bmatrix}
x_{k+1} \\ \lambda_{k}
\end{bmatrix} = \underbrace{
\begin{bmatrix}
A+B\,R^{-1}B^\top A^{-\top}Q & -B\,R^{-1} B^\top A^{-\top} \\
-A^{-\top}Q & A^{-\top}
\end{bmatrix}}_{\mathcal{M}_1}
\begin{bmatrix}
x_{k} \\ \lambda_{k-1}
\end{bmatrix}, \tag{5}
$$
or
$$
\begin{bmatrix}
x_{k+1} \\ \lambda_{k+1}
\end{bmatrix} = \underbrace{
\begin{bmatrix}
A & -B\,R^{-1} B^\top \\
-A^{-\top}Q\,A & A^{-\top}(I + Q\,B\,R^{-1} B^\top)
\end{bmatrix}}_{\mathcal{M}_2}
\begin{bmatrix}
x_{k} \\ \lambda_{k}
\end{bmatrix}. \tag{6}
$$
As $k\to\infty$ both $x_k$ and $u_k$ should go to zero. Therefore $\lambda_k$ should go to zero as well, since $u_k$ is proportional to $\lambda_k$. The states $x_k$ and $\lambda_k$ can only go to zero if the initial conditions only excite stable modes of $\mathcal{M}_1$/$\mathcal{M}_2$, which can be found with the help of the Schur decomposition (the eigen decomposition could be used as well)
$$
\mathcal{M}_i = U_i\,T_i\,U_i^\top, \ \forall\,i\in\{1,2\} \tag{7}
$$
with
$$
U_i =
\begin{bmatrix}
{}^{11}U_i & {}^{12}U_i \\
{}^{21}U_i & {}^{22}U_i
\end{bmatrix}, \tag{8}
$$
$$
T_i =
\begin{bmatrix}
{}^{1}T_i & {}^{2}T_i \\
0 & {}^{3}T_i
\end{bmatrix}, \tag{9}
$$
such that $U_i^\top U_i = I$ and the eigenvalues of ${}^{1}T_i$ are all the eigenvalues of $\mathcal{M}_i$ which have an absolute value smaller than one. In order to be stable the initial conditions should lie in the span of the vectors associated with ${}^{1}T_i$ and thus satisfy
$$
\begin{bmatrix}
x_0 \\ \lambda_{i-2}
\end{bmatrix} =
\begin{bmatrix}
{}^{11}U_i \\ {}^{21}U_i
\end{bmatrix} z_0, \tag{10}
$$
such that they only excite the modes with stable eigenvalues. The value for $x_0$ can not be picked, so in order to still satisfy it one should solve for $\lambda_{i-2}$, which can be shown to yield
\begin{align}
\lambda_{i-2} &= P_i\, x_0, \tag{11} \\
P_i &= {}^{21}U_i\,{}^{11}U_i^{-1}. \tag{12}
\end{align}
By using $(7)$, combined with $U_i^\top U_i = I$ and $(9)$, then the next time step can be shown to yield
$$
\begin{align}
\begin{bmatrix}
x_1 \\ \lambda_{i-1}
\end{bmatrix} &=
U_i\,T_i\,U_i^\top
\begin{bmatrix}
{}^{11}U_i \\ {}^{21}U_i
\end{bmatrix} z_0 \\
&=
\begin{bmatrix}
{}^{11}U_i \\ {}^{21}U_i
\end{bmatrix} {}^{1}T_i\,z_0 \\
\end{align} \tag{13}
$$
From this it can be argued that any future step will remain in the span of $\begin{bmatrix}
{}^{11}U_i^\top & {}^{21}U_i^\top
\end{bmatrix}^\top$ and thus $(11)$ also should hold for all $k>0$.
Substituting $(11)$ into $(4)$ for $i=1$ gives
\begin{align}
u_k &= -R^{-1} B^\top \lambda_k \\
&= -R^{-1} B^\top P_1\, x_{k+1} \\
&= -R^{-1} B^\top P_1\, (A\,x_k + B\,u_k)
\end{align}
Solving this for $u_k$ gives
$$
u_k = -(R + B^\top P_1\,B^\top)^{-1} B^\top P_1\,A\,x_k, \tag{14}
$$
which matches the result from discrete infinite horizon LQR, so $P_1$ should also be the solution to the DARE.
Substituting $(11)$ into $(4)$ for $i=2$ gives
$$
\begin{align}
u_k &= -R^{-1} B^\top \lambda_k \\
&= -R^{-1} B^\top P_2\, x_k
\end{align} \tag{15}
$$
which should still give the same feedback gain as for $i=1$. Calculating it this way might be attractive if $R^{-1}$ can be calculated very quickly compared to $(R + B^\top P_1\,B^\top)^{-1}$, but using this does not directly give a solution to the DARE. The solution could be obtained with for example solve the following discrete Lyapunov equation
$$
(A - B\,K)^\top P\,(A - B\,K) - P + Q + K^\top R\,K = 0, \tag{16}
$$
with $K = R^{-1} B^\top P_2$.
For both $i=1$ and $i=2$ the matrix $\mathcal{M}_i$ uses $A^{-\top}$. This might fail if $A^{-\top}$ is singular. By reformulating is as a generalized eigenvalue problem one can still solve for ${}^{21}U_i$ and ${}^{11}U_i$. Namely the eigenvalue problem $\mu\,\vec{v} = \mathcal{M}_1\,\vec{v}$ associated with $(5)$ can also be written as
$$
\mu \begin{bmatrix}
I & B\,R^{-1} B^\top \\
0 & A^\top
\end{bmatrix}\vec{v} =
\begin{bmatrix}
A & 0 \\
-Q & I
\end{bmatrix}\vec{v}
$$
and the eigenvalue problem $\mu\,\vec{v} = \mathcal{M}_2\,\vec{v}$ associated with $(6)$ as
$$
\mu \begin{bmatrix}
I & -A^\top \\
0 & A^\top
\end{bmatrix}\vec{v} =
\begin{bmatrix}
(Q + I)A & -I-(Q + I)B\,R^{-1} B^\top \\
-Q\,A & I + Q\,B\,R^{-1} B^\top
\end{bmatrix}\vec{v}.
$$
These matrices can then also be used for the generalized Schur decomposition.
Best Answer
You can convert the stated matrix equation into the standard DARE form by using the Woodbury matrix identity, which states
$$ \left(A + U\,C\,V \right)^{-1} = A^{-1} - A^{-1} U \left(C^{-1} + V\,A^{-1} U \right)^{-1} V\,A^{-1}. \tag{1} $$
Using $(1)$ together with $A = I$, $U = F(\gamma) = B_2 B_2^\top - \gamma^{-2} B_1 B_1^\top$, $C = I$ and $V = M$ yields the following
$$ \left(I + F(\gamma)\,M\right)^{-1} = I - F(\gamma)\,(I + M\,F(\gamma))^{-1} M. \tag{2} $$
When assuming that $F(\gamma) \geq 0$ it should also be possible find a $L(\gamma)$ such that $F(\gamma) = L(\gamma)\,L(\gamma)^\top$. One way to find such $L(\gamma)$ would be for example by using the Cholesky decomposition. Substituting this in $F(\gamma)\,(I + M\,F(\gamma))^{-1}$ and partially applying the push-through identity yields
$$ F(\gamma)\,(I + M\,F(\gamma))^{-1} = L(\gamma)\,(I + L(\gamma)^\top M\,L(\gamma))^{-1} L(\gamma)^\top. \tag{3} $$
Substituting this in $(2)$ and its resulting expression in the original matrix equation yields
$$ M = A^\top M\,A - A^\top M\,L(\gamma)\,(I + L(\gamma)^\top M\,L(\gamma))^{-1} L(\gamma)^\top M\,A + C_1^\top C_1, \tag{4} $$
which can be solved for $M$ using a DARE solver using $R = I$ and $B = L(\gamma)$.