Count the number of rational canonical form&find similarity classess
Since you want just to count similarity classes, rather than for instance also find their sizes, and you want to include all matrices, not just invertible ones, it turns out that there is a fairly simple answer to this question, in spite of my pessimistic comment.
You don't really want to classify by the kind of characteristic polynomials, just count the lists of invariant factors one can produce. Remember that one may take the invariant factors to form a finite list of non-constant monic polynomials over $\def\Fq{\Bbb F_q}\Fq$, each one dividing the next, and their product must be of degree$~n$ (it is the characteristic polynomial). One needs to classify these possibilities by the degrees of the polynomials involved. Concretely for a given list of invariant factors let $\mu_k$ be the number of invariant factors of degree at least$~k$, then I will call $\mu=(\mu_1,\mu_2,\ldots)$ the type of the list. It forms a partition of$~n$, with $\mu_1$ the number of invariant factors, and the length$~l(\mu)$ of the partition (the number of nonzero parts) equal to the degree of the largest invariant factor, which is the minimal polynomial. (It might seem more natural to take the degrees of the invariant factors themselves, in reverse (weakly decreasing) order, as type; this is the transpose partition, but the above choice will turn out to have a slight advantage for enumeration).
Now the surprising fact is that the number of possible lists of invariant factors with a given type$~\mu$ depends only on $l(\mu)$ and the order of the field, namely it equals $q^{l(\mu)}$. This is because every invariant factor$~f_i$ must be obtained by taking the previous invariant factor (or the constant $1$ in case of the first invariant factor) and multiplying it by a polynomial with the appropriate degree$~d$; since every monic polynomial of degree is possible, this qives $q^d$ possibilities for choosing$~f_i$. The value of$~d$ is of course the difference of the degree of$~f_i$ and the degree of the previous invariant factor (both are determined by the type), so that the product of the numbers of choices so far becomes $q^{\deg(f_i)}$. In particular, once one has chosen the final invariant factor (the minimal polynomial) the product of the numbers of choices equals $q^{l(\mu)}$. (Note that although this equals the number of monic polynomials of the degree required for the minimal polynomial, it does not mean (unless there is only one invariant factor) that each such polynomial is obtained once as minimal polynomial: some are never obtained (e.g., the irreducible polynomials), while others are obtained as final element of more than one list of invariant factors.)
The number of similarity classes in $M_n(\Fq)$ therefore is $$ \sum_{\mu\in\mathcal P_n}q^{l(\mu)} \qquad\text{where $\mathcal P_n$ is the set of partitions of $n$.} $$ For instance in $M_9(\Fq)$ there are $q +4q^2 +7q^3 +6q^4 +5q^5 +3q^6 +2q^7 +q^8+ q^9$ similarity classes, as can be found by classifying the $30$ partitions of$~9$ by their number of parts; for instance the coefficient $7$ of $q^3$ corresponds to the $7$ partitions $[7,1,1]$, $[6,2,1]$, $[5,3,1]$, $[5,2,2]$, $[4,4,1]$, $[4,3,2]$, $[3,3,3]$ of $9$ into $3$ nonzero parts.
But this result was computed without having to explicitly enumerate partitions, since there is a nice generating series description of these polynomials. If $P_n\in\Bbb Z[X]$ denotes the polynomial such that there are $P_n[q]$ similarity classes in $M_n(\Fq)$, then the coefficient of $X^lY^n$ in the formal series $\sum_{n\in\Bbb N}P_nY^n\in\Bbb Z[X][[Y]]$ counts the number of partitions of $n$ with $l$ nonzero parts, and this series is given by $$ \sum_{n\in\Bbb N}P_nY^n=\prod_{k>0}\frac1{1-XY^k}. $$ For $X:=1$ this gives the well known generating series for the numbers of partitions of$~n$, and since in the expansion of that product choosing a term $Y^{mk}$ in the geometric series $\frac1{1-Y^k}$ corresponds to having $m~$parts equal to$~k$ in the partition, attaching a factor $X^m$ to this term results in recording in the exponent of$~X$ the number of nonzero parts of the partition. So to describe the similarity classes in $M_n(\Fq)$, one may evaluate this product up to the terms of degree$~n$ in$~Y$ (which means only the factors up to $k=n$ of the infinite product need to be considered, and in each the geometric series need only be developed to a limited length), and then take those terms of degree$~n$.