First off, the $2\times2$ case is easy: two distinct $2\times2$ matrices are similar if and only if they have the same characteristic polynomial, and neither of them is a scalar matrix (multiple of the identity). This is essentially because for such small matrices the characteristic polynomial leaves very little room for variation of the minimal polynomial (and other invariant factors that in this case follow from it).

There is a decision method that works in full generality (any field$~K$, any size$~n$ square matrices) and is theoretically simple: compute the invariant factors by determining the Smith normal form of the matrices $XI_n-A$ and $XI_n-B$ (those whose determinants define the respective characteristic polynomials) as matrices with polynomials in$~X$ as entries; then $A$ and $B$ are similar if and only if identical lists of invariant factors are found.

With "theoretically simple" I mean the method only uses basic operations in$~K$ you might be expected to perform: arithmetic of scalars and testing for equality; notably it does not require finding roots or factoring polynomials as any method based on eigenvalues would do. Essentially the Smith normal form algorithm is a generalisation of Gaussian elimination, but working over a PID rather than over a field (which here means: doing polynomial operations on polynomial entries, rather than just using scalars), and allowing column operations as well as row operations. The final goal is to reduce the matrix to its Smith normal form, which has nonzero (polynomial) entries only on the main diagonal, and which are monic polynomials such that each one divides the next (in general there could also be trailing zeros, but that does not happen here since our initial matrices have monic, therefore nonzero, determinants). The diagonal entries found that are different from$~1$ are the invariant factors of $A,B$ respectively; the last one is the minimal polynomial and their product is the characteristic polynomial. Those two polynomials can of course also be computed by other methods, and if either of them differs between $A$ and $B$ this suffices to show they are not similar; in the (rather rare) general case however all invariant factors are needed. In most cases there will be many diagonal entries$~1$, followed by one or very few non-constant polynomials that are the invariant factors.

Being theoretically simple does not mean the procedure is easy to perform (by hand), as it is like in iterated version of the Euclidean algorithm in the ring $K[X]$ of polynomials, which algorithm is notable for producing very complicated intermediate results, at least when $K=\Bbb Q$ (and over $\Bbb R$ or $\Bbb C$ the algorithm isn't really effective, as testing for equality is not accurately possible). It works fine over finite fields though.