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