This is an old prelim problem.
The question is to classify all conjugacy classes of real $4 \times 4$ matrices satisfying $A^3 + A = A^2 + I$. Factoring this gives $(A-I)^2(A+I) = 0$. I can also see that $A(A^2-A +I) = I$, so $A$ is invertible. This question seems to show how to answer questions like this, but I'm not sure I understand the machinery.
It's obvious that if a matrix satisfies this equation, so will any conjugate/similiar matrix. I also understand that we can use $A$ to view $\mathbb{R}^4$ as a $\mathbb{R}[x]$ module where a polynomial $f(x)$ acts on a vector $v$ by $f(x)v = f(A)v$. Since $\mathbb{R}[x]$ is a PID, we can use the classification of finitely generated modules over a PID to write $\mathbb{R}^4$ as a sum of cyclic and free $\mathbb{R}[x]$ modules. I also see how these cyclic factors will give the rational canonical form for the entire conjugacy class, so the problem is really how to find those cyclic factors. If we have all the possible cyclic factors, we can just list all the rational canonical forms and be done.
Is this the right strategy/understanding? How do we find those cyclic factors from this equation, or in general? How does the 4 relate to how many factors we should expect, or is there no relation?
You seem to have the right ideas so far.
In short: a matrix written in its Frobenius normal form will satisfy the equation $p(A) = 0$ if and only if the polynomial associated with each block (i.e. the quotient in each cyclic factor) divides $p(x) = (x - 1)^2 (x + 1)$.
To put it another way: let $C(q)$ denote the companion matrix of the polynomial $q(x)$. For any chain of polynomials $p_1 \mid \cdots \mid p_k \mid p$, the direct sum $$ M = C(p_1) \oplus \cdots \oplus C(p_k) $$ will satisfy $p(M) = 0$. Moreover, every matrix $A$ satisfying $p(A) = 0$ is similar to such a matrix $M$. Note that $M$ will be $4 \times 4$ if and only if the degrees of $p_1,\dots,p_k$ have sum $4$.