An alternative algorithm to find the Jordan form/basis for a complex matrix. I am currently studying System Theory, and the exam involves a lot of finding Jordan forms/bases for state transition matrices. I know there is an algorithm for doing so which involves generalized eigenvectors and all, but that involves computing many powers of the matrix, which is tedious and prone to errors. I seem to have found an alternative method, but I can't find it stated anywhere else and I am wondering why it should be unsound, so I'll drop it here along with an example, and I hope someone will clarify what's wrong with it.
Algorithm
Let's say we have a matrix called $A$, and through the characteristic polynomial we find its eigenvalues $\{\lambda_i\}$ whose algebraic multiplicities are $\{m_i\}$. For each eigenvalue let $T_i$ be $A-\lambda_i I$: through it we can easily find $g_i = \dim \ker T_i$, the geometric multiplicity of the eigenvalue, and the eigenvectors of $A$,  $v_{i,1,1} \cdots v_{i,m_i,1}$, constituting a basis for $\ker T_i$. I'm calling these vectors like this because $v_{i,j,1}$ is the first vector in the partial basis that corresponds to the $j$-th Jordan block with eigenvalue $\lambda_i$.
Then from the shape of the Jordan block we know that $A v_{i,j,k+1} = \lambda_i v_{i,j,k+1} + v_{i,j,k}$ or alternatively $T_i v_{i,j,k+1} = v_{i,j,k}$, so we can recursively find $v_{i,j,k+1}$ through a linear equation, or if we find no solution we know that the Jordan block with partial basis starting with $v_{i,j,1}$ is done.
Then we can put together all the $v_{i,j,k}$'s to find the Jordan form and the desired basis.
Example
We wish to find the Jordan form of
$A=\left(
\begin{array}{cccc}
 2 & 0 & 0 & 0 \\
 0 & 2 & 1 & 0 \\
 0 & 0 & 2 & 0 \\
 1 & 0 & 0 & 2 \\
\end{array}
\right)$.
The characteristic polynomial is $\det(A-\lambda I) = (\lambda - 2)^4$, so we have $\lambda_1 = 2, m_1 = 4$.
We take $\lambda_1 = 2$: by inspecting the kernel of $T_1 = A - 2I$ we find $v_{1,1,1} = \left(
\begin{array}{c}
 0 \\
 1 \\
 0 \\
 0 \\
\end{array}
\right)$ and $v_{1,2,1}=\left(
\begin{array}{c}
 0 \\
 0 \\
 0 \\
 1 \\
\end{array}
\right)$.
We start from $v_{1,1,1}$, and we try to find $v_{1,1,2}$ by solving $T_1 v_{1,1,2} = v_{1,1,1} \implies \left(
\begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & 0 \\
\end{array}
\right) \left(
\begin{array}{c}
 x \\
 y \\
 z \\
 w \\
\end{array}
\right) = \left(
\begin{array}{c}
 0 \\
 1 \\
 0 \\
 0 \\
\end{array}
\right) \implies x=0, z=1$, so I choose $v_{1,1,2}=\left(
\begin{array}{c}
 0 \\
 0 \\
 1 \\
 0 \\
\end{array}
\right)$. Now we try to find $v_{1,1,3}$ by solving $T_1 v_{1,1,3} = v_{1,1,2} \implies \left(
\begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & 0 \\
\end{array}
\right) \left(
\begin{array}{c}
 x \\
 y \\
 z \\
 w \\
\end{array}
\right) = \left(
\begin{array}{c}
 0 \\
 0 \\
 1 \\
 0 \\
\end{array}
\right)$ which is impossible.
So the basis for the first Jordan block is complete, and there must be $v_{1,2,2}$ that is the last vector we seek. We find it as usual by $T_1 v_{1,2,2} = v_{1,2,1} \implies \left(
\begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & 0 \\
\end{array}
\right) \left(
\begin{array}{c}
 x \\
 y \\
 z \\
 w \\
\end{array}
\right) = \left(
\begin{array}{c}
 0 \\
 0 \\
 0 \\
 1 \\
\end{array}
\right) \implies x=1, z=0$ so we choose $v_{1,2,2}=\left(
\begin{array}{c}
 1 \\
 0 \\
 0 \\
 0 \\
\end{array}
\right)$.
So we have found $J=\left(
\begin{array}{cccc}
 2 & 1 & 0 & 0 \\
 0 & 2 & 0 & 0 \\
 0 & 0 & 2 & 1 \\
 0 & 0 & 0 & 2 \\
\end{array}
\right)$ by looking at the lengths of the chains we computed, and also we found $S=\left(
\begin{array}{cccc}
 0 & 0 & 0 & 1 \\
 1 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 \\
 0 & 0 & 1 & 0 \\
\end{array}
\right)$ such that $J=S^{-1}AS$.
This looks long because I wrote it all down, but most times many vectors can be computed mentally, and most importantly it doesn't involve computing $T_1^2$ and $T_1^3$.
 A: I presume that, given the time elapsed since you posed the question, you have already found out the answer. Be it as it may, this is my own comment on your question.
Your algorithm is a good probabilistic (but not deterministic) one. I mean, it may, and probably will, give a right answer to the problem; but it may fail. Let me show this through an example. Suppose you start with a $3 \times 3$ matrix:
$$\begin{pmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix}.$$ Forget that this is already a Jordan matrix. My point is just that the algorithm may fail, so assume that you don't notice and go on with the method.
You compute the characteristic polynomial and find the eigenvalue $1$. The matrix for this eigenvalue will be
$$\begin{pmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0\end{pmatrix}.$$ You compute the space of eigenvectors, which gives the equation $X = 0$, and you choose a basis for this: making a mad election (but my point is just that failure is possible!), you choose the vectors $(0, 1, 2)$ and $(0, 2, 1)$. Then you look for a vector giving $(0, 1, 2)$ by multiplication with the above matrix and find no solution, hence you conclude that this gives one block. You proceed searching for a vector giving $(0, 2, 1)$, but again there is no solution. So, you are left with only two vectors and no basis.
This was just an extreme example. More sophisticated examples could be constructed. The fact is that by arbitrarily choosing the basis for the eigenvectors space you cannot be sure that one particular choice will allow to be extended in the way you want. And with big examples, the failure could appear quite late in the process, forcing you to change the election.
