Solution (source code)

= Solution

Expand $(I+B)^2$ in components:
$$
(A^2)_{ij}=(\delta_{ik}+b_{ik})(\delta_{kj}+b_{kj})=\boxed{\delta_{ij}+2b_{ij}+b_{ik}b_{kj}}.
$$
If $B^2=0$, all powers $B^k$ for $k\geq2$ vanish. The matrix <binomial theorem>, valid here because $I$ commutes with $B$, therefore yields
$$
\boxed{A^n=I+nB,\qquad(A^n)_{ij}=\delta_{ij}+nb_{ij}.}
$$
These <powers of the identity plus a square-zero matrix> can also be proved by induction: $(I+nB)(I+B)=I+(n+1)B$. For the displayed example, $B^2=0$ directly, and
$$
\boxed{\begin{pmatrix}1&p\\0&1\end{pmatrix}^{n}=\begin{pmatrix}1&np\\0&1\end{pmatrix}.}
$$
This is the repeated <function composition> of a <shear mapping>. The argument gives every nonnegative integer $n$; since $(I+B)^{-1}=I-B$, the same formula also holds for negative integers.