[guided]Let $n$ denote the state dimension, let $m$ denote the input dimension, and let $p$ denote the output dimension, so $A,\tilde A,T \in \mathbb{R}^{n \times n}$ with $T$ invertible, $B,\tilde B \in \mathbb{R}^{n \times m}$, and $C,\tilde C \in \mathbb{R}^{p \times n}$. The reachability matrix is built from the columns of $B,AB,\dots,A^{n-1}B$. Thus we need to understand how each block $A^kB$ changes when the state coordinate is replaced by the similar coordinate system determined by $T$.
For $k=0$, the transformed block is just
\begin{align*}
\tilde A^0\tilde B = I_n\tilde B = \tilde B = T^{-1}B = T^{-1}A^0B.
\end{align*}
Now let $k$ be an integer with $1 \le k \le n-1$. Since $\tilde A=T^{-1}AT$, multiplying powers of $\tilde A$ produces adjacent factors $TT^{-1}$, each of which equals $I_n$. Therefore
\begin{align*}
\tilde A^k\tilde B = (T^{-1}AT)^kT^{-1}B.
\end{align*}
Cancelling the adjacent inverse pairs gives
\begin{align*}
(T^{-1}AT)^kT^{-1}B = T^{-1}A^kB.
\end{align*}
So every reachable block in the transformed realization is obtained from the corresponding original block by the same left multiplication by $T^{-1}$. Define the reachability matrices $\mathcal R(A,B) \in \mathbb{R}^{n \times nm}$ and $\mathcal R(\tilde A,\tilde B) \in \mathbb{R}^{n \times nm}$ by
\begin{align*}
\mathcal R(A,B) := [B, AB, \ldots, A^{n-1}B]
\end{align*}
and
\begin{align*}
\mathcal R(\tilde A,\tilde B) := [\tilde B, \tilde A\tilde B, \ldots, \tilde A^{n-1}\tilde B].
\end{align*}
Since the reachability matrix is the horizontal concatenation of these blocks, this gives
\begin{align*}
\mathcal R(\tilde A,\tilde B) = T^{-1}\mathcal R(A,B).
\end{align*}[/guided]