Solve the equation $X^2+X=\text{a given matrix}$ I want to solve the quadratic matrix equation 
$$X^2+X=\begin{pmatrix}1&1\\1&1\end{pmatrix}$$
If I put $X$ in the form 
$$X=\begin{pmatrix}a&b\\c&d\end{pmatrix}$$
then I find complicated equations. Is there a simple way to tackle the problem without using diagonalization?
 A: Write $J=\pmatrix{1&1\\1&1}$. The equation $X^2+X=J$ implies $XJ=JX$. The matrices satisfying $XJ=JX$ are those of the form $X=aI+bJ$. Then $X^2=a^2I+(2ab+2b^2)J$, so we get
$$(a^2+a)I+(b+2ab+2b^2)J=J.$$
Therefore $a^2+a=0$, entailing $a\in\{0,-1\}$ and $b+2ab+2b^2=1$,
which gives a quadratic equation for each of the two possible $a$-values.
A: The spectral mapping theorem gives $$\{\lambda^2 +\lambda : \lambda \in \sigma(X)\} = \sigma(X^2 + X) = \sigma\begin{pmatrix}1&1\\1&1\end{pmatrix} = \{0,2\}$$
It follows that we have one of the four options:


*

*$\sigma(X) = \{0,1\}$ so $\operatorname{Tr} X = 1$.


A calculation with $X = \pmatrix{a & b \\ c & 1-a}$ yields
$$\pmatrix{a+a^2+bc & 2b \\ 2c & a^2-3a+bc+2} = X^2 + X = \begin{pmatrix}1&1\\1&1\end{pmatrix}$$
so the only solution is $a = b = c = \frac12$.


*

*$\sigma(X) = \{0,-2\}$ so $\operatorname{Tr} X = -2$.


A calculation with $X = \pmatrix{a & b \\ c & -2-a}$ yields
$$\pmatrix{a+a^2+bc & -b \\ -c & a^2+3a+bc + 2} = X^2 + X = \begin{pmatrix}1&1\\1&1\end{pmatrix}$$
so the only solution is $a = b = c = -1$.


*

*$\sigma(X) = \{-1,-2\}$ so $\operatorname{Tr} X = -3$.


A calculation with $X = \pmatrix{a & b \\ c & -3-a}$ yields
$$\pmatrix{a+a^2+bc & -2b \\ -2c & a^2+5a + bc + 6} = X^2 + X = \begin{pmatrix}1&1\\1&1\end{pmatrix}$$
so $a = -\frac32$, $b = c = -\frac12$ is the only solution.


*

*$\sigma(X) = \{-1,1\}$ so $\operatorname{Tr} X = 0$.


A calculation with $X = \pmatrix{a & b \\ c & -a}$ yields
$$\pmatrix{a+a^2+bc & b \\ c & a^2-a+bc} = X^2 + X = \begin{pmatrix}1&1\\1&1\end{pmatrix}$$
so $a = 0$, $b = c = 1$ is the only solution.
So the solutions for $X$ are
$$\pmatrix{\frac12 & \frac12 \\ \frac12 & \frac12},\quad \pmatrix{-1 & -1 \\ -1 & -1},\quad \pmatrix{-\frac32 & -\frac12 \\ -\frac12 & -\frac32},\quad \pmatrix{0 & 1 \\ 1 & 0}$$
A: Your method also works, though cumbersome (not asap):
$$X^2+X=1 \Rightarrow \begin{pmatrix}a&b\\c&d\end{pmatrix}^2+\begin{pmatrix}a&b\\c&d\end{pmatrix}=\begin{pmatrix}1&1\\1&1\end{pmatrix} \Rightarrow\begin{cases}a^2+a+bc=1\\ ab+bd+b=1\\ ac+c+cd=1\\ bc+d^2+d=1\end{cases}.$$
Subtract the last from the first:
$$a^2-d^2+a-d=0 \Rightarrow (a-d)(\underbrace{a+d+1}_{\ne 0 \ \text{from (2)}})=0 \Rightarrow d=a.$$
Plug $d=a$ in the second and third equations and subtract:
$$(2ab+b)-(2ac+c)=0 \Rightarrow (b-c)(\underbrace{2a+1}_{\ne 0 \ \text{from (2)}})=0 \Rightarrow c=b.$$
Plug $c=b$ into the first and second equations and add them: 
$$a^2+a+b^2+2ab+b=2 \Rightarrow (a+b)^2+(a+b)-2=0 \Rightarrow \\
a+b=-2;1 \Rightarrow 1) \ b=-a-2; \ 2) \ b=-a+1.$$
$1$) Plug $b=-a-2$ into the first equation:
$$a^2+a+(-a-2)^2=1 \Rightarrow 2a^2+5a+3=0 \Rightarrow a=-\frac32; -1 \Rightarrow \\
X_1=\begin{pmatrix}a&b\\c&d\end{pmatrix}=\begin{pmatrix}-\frac32&-\frac12\\-\frac12&-\frac32\end{pmatrix};\\
X_2=\begin{pmatrix}a&b\\c&d\end{pmatrix}=\begin{pmatrix}-1&-1\\-1&-1\end{pmatrix}.$$
$2$) Plug $b=-a+1$ into the first equation:
$$a^2+a+(-a+1)^2=1 \Rightarrow 2a^2-a=0 \Rightarrow a=0; \frac12 \Rightarrow \\
X_3=\begin{pmatrix}a&b\\c&d\end{pmatrix}=\begin{pmatrix}0&1\\1&0\end{pmatrix};\\
X_4=\begin{pmatrix}a&b\\c&d\end{pmatrix}=\begin{pmatrix}\frac12&\frac12\\\frac12&\frac12\end{pmatrix}.$$
A: The following SymPy script
from sympy import *

x11, x12, x21, x22 = symbols('x11 x12 x21 x22', real=True)

X = Matrix([[ x11, x12],
            [ x21, x22]])

R = X**2 + X - ones(2,2)

solutions = solve_poly_system([R[0,0], R[0,1], R[1,0], R[1,1]], x11, x12, x21, x22)

for (a,b,c,d) in solutions:
    print Matrix([[a,b],
                  [c,d]])

outputs the following $4$ solutions
Matrix([[-3/2, -1/2],
        [-1/2, -3/2]])
Matrix([[-1, -1],
        [-1, -1]])
Matrix([[0, 1],
        [1, 0]])
Matrix([[1/2, 1/2],
        [1/2, 1/2]])

A: Using Jordan canonical form, $X$ must be similar to one of the three matrices
\begin{align}
\begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
\end{pmatrix},
\begin{pmatrix}
\lambda & 0 \\
0 & \lambda
\end{pmatrix},
\begin{pmatrix}
\lambda & 1 \\
0 & \lambda
\end{pmatrix}
\end{align}
, where $\lambda_1 \neq \lambda_2$, which implies that $X^2 + X$ must be similar to 
\begin{align}
\begin{pmatrix}
\lambda_1^2 + \lambda_1 & 0 \\
0 & \lambda_2^2 + \lambda_2
\end{pmatrix},
\begin{pmatrix}
\lambda^2 + \lambda & 0 \\
0 & \lambda^2 + \lambda
\end{pmatrix},
\begin{pmatrix}
\lambda^2 + \lambda & 2\lambda + 1 \\
0 & \lambda^2 + \lambda
\end{pmatrix}
\end{align}
, respectively. Since the latter two matrices have only a single eigenvalue $\lambda^2 + \lambda$, it is impossible for them to be similar to $J \equiv \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$, which has two distinct eigenvalues. Therefore, $X$ must take the form 
$$X = P^{-1}\begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
\end{pmatrix}P, \tag{1}$$
where $P := \begin{pmatrix} a & b \\ c & d \end{pmatrix}$ is invertible. In terms of $(1)$, the equation $X^2 + X = J$ is expressed as 
$$P^{-1}\begin{pmatrix}
\lambda_1^2 + \lambda_1 & 0 \\
0 & \lambda_2^2 + \lambda_2
\end{pmatrix}P = J. \tag{2}$$
Because $J$ has two distinct eigenvalues $0$ and $2$, without of loss of generality, we can set 
$$\begin{cases} 
\lambda_1^2 + \lambda_1 = 0 \\
\lambda_2^2 + \lambda_2 = 2,
\end{cases} \tag{3}$$
under which $(2)$ is equivalent to 
$$\begin{pmatrix}
0 & 0 \\
0 & 2
\end{pmatrix}
\begin{pmatrix} a & b \\ c & d \end{pmatrix} = 
\begin{pmatrix} a & b \\ c & d \end{pmatrix}\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}, $$
or
$$\begin{pmatrix}
0 & 0 \\
2c & 2d
\end{pmatrix} = 
\begin{pmatrix}
a + b & a + b \\
c + d & c + d
\end{pmatrix}, $$
which requires $a = -b$ and $c = d$ (hence $a \neq 0, c \neq 0$). Therefore, 
$$P = \begin{pmatrix}
a & -a \\
c & c 
\end{pmatrix}, \quad
P^{-1} = \frac{1}{2ac}\begin{pmatrix}
c & a \\
-c & a
\end{pmatrix}. \tag{4}$$
Substituting $(4)$ into $(1)$ gives
$$X = \frac{1}{2}\begin{pmatrix}
\lambda_1 + \lambda_2 & \lambda_2 - \lambda_1 \\
\lambda_2 - \lambda_1 & \lambda_1 + \lambda_2
\end{pmatrix}. \tag{5}$$
While $(3)$ admits the following four combinations of $(\lambda_1, \lambda_2)$:
$$(0, -2), (0, 1), (-1, -2), (-1, 1),$$
$(5)$ entails the following four solutions of $X$:
\begin{align}
\begin{pmatrix}
-1 & -1 \\
-1 & -1
\end{pmatrix},
\begin{pmatrix}
1/2 & 1/2 \\
1/2 & 1/2
\end{pmatrix},
\begin{pmatrix}
-3/2 & -1/2 \\
-1/2 & -3/2
\end{pmatrix},
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}.
\end{align}
A: There is a matrix $S$ such that
$$
J=S\begin{bmatrix} 2 & 0 \\ 0 & 0 \end{bmatrix}S^{-1}
$$
If you set $Y=S^{-1}XS$, then $X=SYS^{-1}$ and the equation becomes
$$
S(Y^2+Y)S^{-1}=S\begin{bmatrix} 2 & 0 \\ 0 & 0 \end{bmatrix}S^{-1}
$$
that simplifies to
$$
Y^2+Y=\begin{bmatrix} 2 & 0 \\ 0 & 0 \end{bmatrix}
$$
We can try and complete the square:
$$
4Y^2+4Y=\begin{bmatrix} 8 & 0 \\ 0 & 0 \end{bmatrix}
$$
so the equation becomes
$$
(2Y+I)^2=\begin{bmatrix} 9 & 0 \\ 0 & 1 \end{bmatrix}
$$
Now the problem is simpler; if $2Y+I=\left[\begin{smallmatrix}a&b\\c&d\end{smallmatrix}\right]$, then the equation becomes
$$
\begin{bmatrix}
a^2 + bc & ab + bd \\
ac + cd & bc + d^2
\end{bmatrix}=
\begin{bmatrix} 9 & 0 \\ 0 & 1 \end{bmatrix}
$$
Case 1: $b=0$; then $a=\pm3$, $d=\pm1$ and $c=0$
Case 2: $a+d=0$: then $a^2+bc=9$ and $a^2+bc=1$: a contradiction.
Thus we have four solutions and
$$
Y=\frac{1}{2}\left(\begin{bmatrix} \pm3 & 0 \\ 0 & \pm1 \end{bmatrix}-I\right)
$$
so
$$
X=\frac{1}{2}S\left(\begin{bmatrix} \pm3 & 0 \\ 0 & \pm1 \end{bmatrix}-I\right)S^{-1}
$$
How do you determine the matrix $S$? Just determine an eigenvector relative to $2$ and one relative to $0$:
$$
S=\begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}
$$
A: For those amongst you folks, like egreg and myself, who enjoy completing the square, we can also have at it like this:
Set
$J = \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}; \tag 1$
then our equation reads
$X^2 + X = J, \tag 2$
from which we may write, completing the square, 
$(X + \dfrac{1}{2}I)^2 = X^2 + X + \dfrac{1}{4} = J + \dfrac{1}{4}I = \begin{bmatrix} \dfrac{5}{4} & 1 \\ 1 & \dfrac{5}{4} \end{bmatrix}; \tag 3$
our task is now to find the matrices $C$ such that
$C^2 = J + \dfrac{1}{4}I = \begin{bmatrix} \dfrac{5}{4} & 1 \\ 1 & \dfrac{5}{4} \end{bmatrix}; \tag 4$
to this end we observe that, setting
$P = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \tag 5$
whence
$P^2 = I, \tag 6$
we have
$\left (\dfrac{1}{2}I + P \right )^2 = \dfrac{1}{4} I^2 + 2\left ( \dfrac{1}{2}I \right ) P + P^2 = \dfrac{1}{4} I + I + P =  \begin{bmatrix} \dfrac{5}{4} & 1 \\ 1 & \dfrac{5}{4} \end{bmatrix}, \tag 7$
and also
$\left (I + \dfrac{1}{2} P \right )^2 = I^2 + 2I \left (\dfrac{1}{2}P \right ) + \dfrac{1}{4} P^2 = I + \dfrac{1}{4}I + P =  \begin{bmatrix} \dfrac{5}{4} & 1 \\ 1 & \dfrac{5}{4} \end{bmatrix}; \tag 8$
we are thus motivated to generalize by seeking solutions of the form $aI + bP$:
$(aI + bP)^2 = a^2 I^2 + 2abIP + b^2 P^2 = a^2 I + b^2 I + 2ab P$
$= (a^2 + b^2)I + 2ab P =  \begin{bmatrix} \dfrac{5}{4} & 1 \\ 1 & \dfrac{5}{4} \end{bmatrix} = \dfrac{5}{4}I + P; \tag 9$
thus,
$a^2 + b^2 = \dfrac{5}{4}; \; 2ab = 1; \tag{10}$
it follows from (10) that the ordered pairs $(a, b)$ lie on both the circle of radius $5/4$ centered at the origin in the $(a, b)$ plane, and also on the hypberbola $ab = 1/2$; since the points $\pm (1/\sqrt 2, 1/\sqrt 2)$ lie in the interior of the circle, the branches of this hyperbola cross the circle at exactly four points, and we have seen in (8)-(9) that two of them are $(1, 1/2)$ and $(1/2, 1)$; it is then evident from the symmetry $(a, b) \leftrightarrow (-a, -b)$ of the solutions to (10) that there are a total of precisely four pairs $(a, b)$ which solve (10) and hence (9), viz.
$(a, b) = \pm \left (1, \dfrac{1}{2} \right ), \; \pm \left ( \dfrac{1}{2}, 1 \right ); \tag{11}$
since
$X + \dfrac{1}{2}I = aI + bP = \pm \begin{bmatrix} 1 & \dfrac{1}{2} \\ \dfrac{1}{2} & 1 \end{bmatrix}, \pm \begin{bmatrix} \dfrac{1}{2} & 1 \\ 1 & \dfrac{1}{2} \end{bmatrix}, \tag{12}$
we have
$X = \begin{bmatrix}  \dfrac{1}{2} & \dfrac{1}{2} \\ \dfrac{1}{2} &  \dfrac{1}{2} \end{bmatrix}, \; \begin{bmatrix} -1 & -1 \\ -1 & -1  \end{bmatrix}, \; \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \begin{bmatrix} -\dfrac{3}{2} & -\dfrac{1}{2} \\ -\dfrac{1}{2} & -\dfrac{3}{2} \end{bmatrix}, \tag{13}$
as the solutions to (2).
One question does, however, remain:  how do we know every $C$ as in (4) may be cast in the form $aI + bP$?  Well, as pointed out by Lord Shark the Unknown in his answer, $XJ = JX$ implies
$X = aI + bJ; \tag{14}$
but
$J = I + P, \tag{15}$
so in fact
$a I + b J = aI + b(I + P) = (a + b)I + bP, \tag{16}$
which shows that $X$ may be represented in the form $aI + bP$ as well; thus so may
$C = X + \dfrac{1}{2}I = (a + \dfrac{1}{2})I + bJ = (a + b + \dfrac{1}{2})I + bP. \tag{17}$
