Compute a matrix $P$ such that $P$ satisfies $3I+P+P^2=\left(\begin{smallmatrix}3&0&0\\3&6&0\\0&0&6\end{smallmatrix}\right)$ Question: Compute a matrix $P$ such that the matrix $P$ satisfies $3I + P + P^2 = \begin{pmatrix}  3 & 0 & 0\\ 3 & 6 & 0 \\ 0 & 0 & 6 \end{pmatrix}$
Attempt: First, I simplify the equation and get $P + P^2 = P(I + P) = \begin{pmatrix} 0 & 0 & 0 \\3 & 3 & 0\\ 0& 0 & 3\end{pmatrix}$. And then I don't know how to solve further. I think this question somehow is related to the characteristic polynomial, but I am not sure. I get stuck on this question for hours, could someone please help me out, thanks!
 A: This is a bit of a brute-force diagonalisation approach - I do feel as though there should be a simpler method, but this at least gives you an answer.

We can reduce this to solving an equation of the form $X^2=C$ by completing the square, as we would do in the real case. Note that (dropping the now-implicit identity matrix for scalars)
$$\left(P + \frac{1}{2}\right)^2 = P^2 + P + \frac{1}{4} \implies P^2 + P = \left(P + \frac{1}{2}\right)^2 - \frac{1}{4}$$
Writing $Q$ for the RHS of your original equation, we have
$$P^2 + P + 3 = Q \iff \left(P + \frac{1}{2}\right)^2 = Q - \frac{11}{4}$$
We now seek a matrix $X$ such that $X^2 = Q'$ where $Q' = Q - 11/4$. One way to do this is to diagonalise $Q'$ as
$$Q' = V^{-1} D V$$
where $D$ is diagonal, since then we find that taking
$$X = V^{-1} \sqrt{D} V \implies X^2 = Q'$$
and $\sqrt{D}$ is easy to calculate since we just take the square root of the diagonal entries. Indeed, such a decomposition is possible, with
$$V = \begin{pmatrix} -1 & 0 & 0 \\ 0 & 0 & 1 \\ 1 & 1 & 0 \end{pmatrix}, \quad D = \frac{1}{4} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 13 & 0 \\ 0 & 0 & 13 \end{pmatrix}$$
Computing $X$, as indicated above, one can show that
$$
X = \frac{1}{2} \begin{pmatrix} 1 & 0 & 0 \\ \sqrt{13}-1 & \sqrt{13} & 0 \\ 0 & 0 & \sqrt{13} \end{pmatrix}$$
and hence
$$P = \frac{1}{2} \begin{pmatrix} 0 & 0 & 0 \\ \sqrt{13}-1 & \sqrt{13}-1 & 0 \\ 0 & 0 & \sqrt{13}-1 \end{pmatrix} = \frac{\sqrt{13}-1}{2} \begin{pmatrix} 0 & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}
$$
A: Let $C$ be the matrix on the RHS and $f(x)=x^2+x+3$. We want to solve $f(P)=C$.
Since $C$ is the direct sum of two diagonalisable submatrices $\pmatrix{3&0\\ 3&6}$ and $6$, it can be diagonalised as $C=MDM^{-1}$ where $D=\operatorname{diag}(3,6,6)$. (There is no need to calculate $M$; see below.)
Pick two roots $a$ and $b$ of $f(a)=3$ and $f(b)=6$ respectively. Let $g$ be the linear polynomial such that $g(3)=a$ and $g(6)=b$. One may obtain $g$ from the point-slope form $g(x)-a=\frac{b-a}{6-3}(x-3)$.
Then $f(g(D))=f(\operatorname{diag}(a,b,b))=D$. Therefore
$$
f(g(C))=f(g(MDM^{-1}))=Mf(g(D))M^{-1}=MDM^{-1}=C,
$$
i.e. $P=g(C)$ is a solution to $f(P)=C$.
