Euler's Theorem, but stronger

The max order of a matrix in $GL(n,p)$ is $a_n = p^n-1.$ The sum of these is $(p^{p+1} + p^p + \cdots + p^3 + p^2 + p) - (p+1) = (p^{p+1} + p^p + \cdots + p^3) + p(p-1)+1 \cdot (p-1),$ so in base $p$ this is $11\cdots10xx,$ where $x=p-1.$ There are $p-1$ 1's, so the sum is $3(p-1) = \fbox{6048}.$

Perhaps you are looking for a proof that the max order is $p^n-1$ ? Well, first of all, there is a matrix whose order is $p^n-1$ : look at the finite field ${\mathbb F}_{p^n}$ as a vector space of dimension $n$ over ${\mathbb F}_p,$ and look at the matrix of multiplication by a primitive root : its order is $p^n-1$ by definition. On the other hand, the max order is at most $p^n-1,$ because every power of $A$ can be represented as a linear combination of $I,A,A^2,\ldots,A^{n-1}$ by Cayley-Hamilton , and there are exactly $p^n-1$ nonzero such linear combinations, so two of the elements of the list $I,A,A^2,\ldots,A^{p^n-1}$ must coincide.