Show that the sequence $X_n$ converges to a limit $Y$ Let $A$ be a positive definite $n\times n$ matrix. We use the iteration of the mapping $$f(X)=0.5(X^2+B), \ X\in \mathbb{R}^{n\times n}$$ where $B=I-A$.
Show that the sequence $X_0:=0$ (zero matrix), $X_1=f(X_0), X_2=f(f(X_0)), \ldots $ converges to a limit $Y$ if $\|B\|<1$ and show that $(I-Y)(I-Y)=A$.
We want to show that for every $\epsilon>0$ there is a natural number $N$ such that $|X_n-Y|<\epsilon$, right? How could we do that?
 A: We have $$f(X) = \frac 12X^2 + \frac 12 B$$
Let's try to see what initial terms of $(X_n)$ look like (with $X_0=0$):
$$\begin{align}X_1 &= f(0) = \frac 12 0^2 + \frac 12 B = \frac 12 B\\[2mm]
X_2 &= f(f(0)) = \frac12 \bigg(\frac 12 B\bigg)^2 + \frac 12 B = \frac{1}{2^3}B^2 + \frac{1}{2}B \\[2mm]
X_3 &=f(f(f(0)))= \frac 12\bigg( \frac{1}{2^3}B^2 + \frac{1}{2}B \bigg)^2 + \frac 12 B = \frac1{2^7}B^4 + \frac{1}{2^5}B^3 + \frac {1}{2^3}B^2 + \frac{1}{2} B\end{align}$$
We now notice that (which can be proven by induction):
$$X_n = \frac{1}{2} B + \frac {1}{2^3}B^2 + \dots + \frac{1}{2^{2a_n - 1}}B^{a_n} = \sum_{k=1}^{a_n}\frac{1}{2^{2k - 1}}B^{k}$$
where $a_n = a_{n-1}^\text{th}~ \text{triangular number} + 1$, or explicitly
$$\quad a_{n+1} = \frac{a_n(a_n+1)}{2} + 1$$
for all $n \in \mathbb N$ and $a_0 = 1$.

Why exactly this number? This is because the expression $(c_1 + \dots + c_n)^2$ has $\frac{n(n+1)}{2}$ (which is $n^\text{th}$ triangular number) summands, and we add $1$ to  this number since the function $f$ is adding one extra "constant" term.

Now, since $\|B\| < 1$, the following series
$$\lim_{n\to\infty}X_n = \sum_{n=0}^\infty \frac{1}{2^{2k - 1}}B^{k}$$
is convergent.

Finally, let the limit be $Y$. So, we have
$$X_{n+1} = f(X_n)$$
Now "limiting" both sides, we get
$$Y = f(Y) \\
Y = \frac 12 Y^2 + \frac 12 B \\
Y^2 - 2Y + B = 0$$
By the original substitution $B = I-A$:
$$Y^2 - 2Y + I = A \\
(Y - I)^2 = A \\
(I-Y)(I-Y) = A$$
as desired.
A: By mathematical induction, we see that each $X_k$ is a polynomial in $A$. Since $A$ is diagonalisable, it suffices to prove the problem statement in the scalar case.
Presumably, the matrix norm in question is submultiplicative. Hence $\rho(B)\le\|B\|<1$ and in the scalar case, we may assume that $|B|<1$. From $X_0\le|B|<1$, we can prove by mathematical induction that $X_k=\frac12(X_{k-1}^2+B)\le|B|$. Clearly we also have $\frac12(X_{k-1}^2+B)\ge\frac12 B\ge-|B|$. Therefore $|X_k|\le|B|$ for every $k$ and
$$
|X_{k+1}-X_k|
=\left|\frac{X_k^2-X_{k-1}^2}{2}\right|
=\left|\frac{X_k+X_{k-1}}{2}\right|\left|X_k-X_{k-1}\right|
\le|B|\left|X_k-X_{k-1}\right|.
$$
Hence $\{X_k\}$ is a Cauchy sequence and it converges to some $Y\in\mathbb R$. By passing both sides of $X_k=\frac12(X_{k-1}^2+B)$ to the limit, we get $Y=\frac12(Y^2+B)$, which can be rewritten as $(1-Y)^2=A$.
