In order to find the order of a matrix over Z, you can look at the eigenvalues of the matrix.. If the matrix isn't diagonalizable, or if it has an eigenvalue that is not a root of unity, then its order is infinite. Otherwise, the order of the matrix is the LCM of the orders of the roots of unity.
Note that GL(3,Z) has infinitely many elements of order 6, because if A is any matrix of order 6 and B=[[1,1,0],[0,1,0],[0,0,1]] then {B^-nAB^n: n in Z} is an infinite set of matrices of the same order.
For finite groups you could use Gap's conjugacy class functionality. This is not directly exposed in Sage, so you end up talking "Gap". You might be better off asking Gap directly, but it's not too bad (and if you start out in Sage, you get the benefit that you can let Sage do the conversions for you)
sage: G=GL(3,GF(5))
sage: gapG=G._gap_()
sage: CG=gapG.ConjugacyClasses()
sage: sum(c.Size() for c in cc if c.Representative().Order() == 6)
46500
I suspect there are much better representation-theoretic methods of arriving at this number.