Note: for the notation I'm used to, the 1s for $C$ are subdiagonal, as in Wikipedia, not superdiagonal, as in your question.
Under the assumption that $A$ is conjugate to a companion matrix:
If you are willing to accept a probabilistic answer, then there is a very efficient algorithm. Choose $\vec{v} \in \mathbb{R}^n$ randomly, under some reasonable random distribution (mainly: that it doesn't have its support only in a subvariety of $\mathbb{R}^n$). Then let $P$ be the array that consists of $\vec{v}, A\vec{v}, A^2 \vec{v}, ..., A^{n - 1} \vec{v}$. Let $Q = P^{-1}$. Then I claim that $A = Q^{-1} C Q$.
The above claim is equivalent to the following claim: that almost every vector $\vec{v}$ is cyclic, that is, that $\{A^i \vec{v}\}_{i = 0}^{n - 1}$ is linearly independent (formally: that the set of noncyclic vectors is a subvariety of codimension 1).
Proof that $A = Q^{-1} C Q$:
If $\{A^i \vec{v}\}_{i = 0}^{n - 1}$ is linearly independent, then it is a basis of $\mathbb{R}^n$. Therefore, we only need to prove that $A (A^i \vec{v}) = Q^{-1} C Q A^i \vec{v}$ for $0 \leq i \leq n - 1$. But by the definition of $Q$, we have that $Q A^i \vec{v} = \vec{e}_i$, the $i$th standard basis element. For $i \neq n - 1$, we have that $C \vec{e}_i = \vec{e}_{i + 1}$; for $i = n - 1$, w have $C \vec{e}_{n - 1} = \sum_{j = 0}^{n - 1} -u_j \vec{e}_j$.Then for $i \neq n - 1$, we have:
$Q^{-1} C Q A^i \vec{v} = Q^{-1} C \vec{e}_i = Q^{-1} \vec{e}_{i + 1} = A^{i + 1} \vec{v} = A A^i \vec{v}$
and for $i = n - 1$, we have
$Q^{-1} C Q A^{n - 1} \vec{v} = Q^{-1} C \vec{e}_{n - 1} = Q^{-1} \sum_{j = 0}^{n - 1} -u_j \vec{e}_j = \sum_{j = 0}^{n - 1} -u_j A^j \vec{v}$. But by the fact that $\chi_A(x) = x^n + \sum_{j = 0}^{n - 1} u_j x^j$, and the fact that $\chi_A(A) = 0$, we have that $\sum_{j = 0}^{n - 1} -u_j A^j = A^n$. Therefore, $Q^{-1} C Q A^{n - 1} \vec{v} = A^n \vec{v} = A A^{n - 1} \vec{v}$.
Therefore, we've shown that all the basis elements go where they should; therefore, we have that $A = Q^{-1} C Q$.
Proof that almost every vector is cyclic: The assumption that $A$ is conjugate to a companion matrix is equivalent to the assumption that there is some cyclic vector, as the companion matrix has cyclic vector $\vec{e}_1$. The set of non-cyclic vectors is determined by the equation $det(P_\vec{v}) = 0$. As there is some cyclic vector, we have that $det(P_\vec{v})$ is nonzero somewhere, so the set where it is 0 is not the entire space, and therefore is a subvariety of codimension 1.
Corrected according to comment: This determines such a $Q$ in $O(n^3)$ time as applying a matrix to a vector takes $O(n^2)$ and we do this $n$ times, then invert, which takes less than $O(n^3)$. This is as fast as the corresponding method for Jordan normal form according to What is the best algorithm to find the smallest nonzero Eigenvalue of a symmetric matrix? . It has the advantage of working in any field without having to go through any field extensions - which you seemed to indicate was a problem when talking about complex numbers.
Best Answer
I don't quite follow the first paragraph of your question. But reading the rest of your post, if you only need a set $S$ of $N$ $m$-dimensional vectors over the finite field of some small order in which any subset $S' \subset S$ of cardinality $m$ is a set of linearly independent vectors, that's a parity-check matrix of an MDS code over $\mathbb{F}_q$.
Take an $m \times N$ parity-check matrix $H$ of an $[N, N-m, m+1]_q$ code (i.e., a $q$-nary linear code of length $N$, dimension $N-m$, and minimum distance $m+1$). Then whichever $m$ or fewer columns you pick from $H$, they're always linearly independent.
The binary case is no good because the only MDS codes are the trivial ones. For larger $q$, there are known nontrivial MDS codes. Reed-Solomon codes are good examples. They're cyclic codes so you can realize the codewords as an ideal of the polynomial ring $\mathbb{F}_q[x]/(x^N-1)$; you can generate them through multiplication between a certain monic polynomial that divides $x^N-1$ and all polynomials of degree less than $N-m$ (including $0$). Maybe this is systematic enough to work for your purpose?
In any case, you might want to clarify your question a bit or, if MDS codes are indeed what you would like to construct, pick your favorite coding theory textbook and read up on them.