Interesting matrix construction question 
We are told that an $n \times n$ matrix $P$ satisfies $P^3=P$. Can we construct such a matrix $P$?

Of course, we see that $-1,0,1$ are its only eigenvalues and, thus, a diagonal matrix with diagonal elements $0,1,-1$ would do. However, if I were told that the given matrix should not have $0,1,-1$ as its diagonal elements, then what would be the way out?
 A: Any matrix $P$ satisfying $P^3 = P$ satisfies $f(P) = 0$, where $f(x) = (x + 1) x (x - 1)$, so the minimal polynomial $m_p$ of $P$ divides $(x + 1) x (x - 1)$. In particular, $m_p(x)$ has distinct linear factors, equivalently, it is diagonalizable. Now, appeal to the Jordan normal form.

 For any solution $P$ of $P^3 = P$ there is some invertible matrix $Q$ satisfying $$P = Q \pmatrix{\lambda_1\\&\ddots\\&&\lambda_n} Q^{-1},$$ where (per the observation in the question, or using our observation that $m_p \mid (x + 1) x (x - 1)$)  $\lambda_i \in \{-1, 0, 1\}$, $i = 1, \ldots, n$. Conversely, substituting shows that any $P$ of this form satisfies $P^3 = P$, and so matrices of that form comprise precisely the solutions of the equation.

A: The point of this answer is to show that such a $P$ is diagonalisable.  
Let $P$ represent a linear transformation such that $P^{n+1}=P$. If we denote the $n$-th roots of unity by$$\omega_{n,k}=\exp\left(2\pi i\frac{k}{n}\right)$$
then obviously the only eigenvalues of $P$ can be the $\omega_{n,k}$ and $0$.  
We can write the identity transformation as as:
$$
\mathrm{Id}=\mathrm{Id}-P^n+\sum_{r=0}^{n-1}R_r
$$
where
$$
R_r=\frac{1}{n}\sum_{k=0}^{n-1}\omega_{n,k}^rP^{n-k}
$$
Indeed:
$$
\sum_{r=0}^{n-1}R_r=\sum_{r=0}^{n-1}\frac{1}{n}\sum_{k=0}^{n-1}\omega_{n,k}^rP^k=\sum_{k=0}^{n-1}\left(\frac{1}{n}\sum_{r=0}^{n-1}\omega_{n,k}^r\right)P^{n-k}=\sum_{k=0}^{n-1}\delta_{k,0}P^{n-k}=P^n
$$
For the second equality I summed the geometric series when $k\neq0$, while $\delta$ is the Kronecker delta.  
At the same time, notice that
\begin{align}
PR_r&=\frac{1}{n}\sum_{k=0}^{n-1}\omega_{n,k}^rP^{n-k+1}=\frac{1}{n}\omega_{n,r}\sum_{k=1}^{n-1}\omega_{n,k-1}^rP^{n-(k-1)}+\frac{1}{n}\omega_{n,0}^r
P^{n+1}=\\
&=\frac{1}{n}\omega_{n,r}\sum_{k=0}^{n-2}\omega_{n,k}^rP^{n-k}+\frac{1}{n}\omega_{n,r}\,\omega_{n,n-1}^rP^{n-(n-1)}=\omega_{n,r}\frac{1}{n}\sum_{k=0}^{n-1}\omega_{n,k}^rP^{n-k}=\omega_{n,r}R_r
\end{align}
And:
$$
R_r^2=\frac{1}{n^2}\sum_{k=0}^{n-1}\sum_{l=0}^{n-1}(\omega_{n,k}\omega_{n,l})^rP^{2n-k-l}=\frac{1}{n^2}\sum_{k=0}^{n-1}\sum_{l=0}^{n-1}\omega_{n,k+l}^rP^{n-(k+l-n)}=\frac{1}{n}\sum_{m=0}^{n-1}\omega_{n,m}^rP^{n-m}=R_r
$$
Here I used the fact that there are exactly $n$ ways to choose a pair $(k,l)$ such that $k+l=x\,\,\mathrm{mod}\,\,n$ for any given $0\leq x\leq n-1$ and also the facts that $\omega_{n,k+l}=\omega_{n,(k+l\,\,\mathrm{mod}\,\,n)}$ and $P^{n-(k+l-n)}=P^{n-(k+l\,\,\mathrm{mod}\,\,n)}$ (all this is perhaps a bit finicky to show, but it does work - if you want to convince yourself, the simplest way is probably to treat seperately the cases when $k+l\geq n$ and $k+l< n$).
Also, we have $P(\mathrm{Id}-P^n)=0$. If we take some $v$, we find:
$$
v=(\mathrm{Id}-P^n)v+\sum_{r=0}^{n-1}R_rv
$$
This is a decomposition of an arbitrary vector into eigenvectors (the $(\mathrm{Id}-P^n)v$ and $R_rv$), and so if we pick arbitrary bases for the eigenspaces $E_P(0)$ and $E_P(\omega_{n,r})$, these together form a basis for the entire space. Therefore, $P$ is diagonalisable and $\mathrm{Id}-P^n,R_r$ are its associated projectors. Now you can proceed as Giuseppe Negro suggested, by conjugating appropriate diagonal matrices with invertible ones. This should give all possible matrices such that $P^{n+1}=P$.
