To find the geometric multiplicity of an eigenvalue $\lambda$, you want to find
$nullity(A-\lambda I)=n-rank(A-\lambda I)$, where n is the number of columns of A.
(Notice that you want to use n instead of rank(A).)
Any basis of $\ker(A-\lambda I)$ will give you a set of linearly independent eigenvectors for the eigenvalue $\lambda$.
In your example, you can find a generalized eigenvector w for $\lambda=2$ by either selecting an eigenvector v and then solving $(A-2I)w=v$ for w, or by choosing any vector w which is not in $\ker(A-2I)$ and then taking $v=(A-2I)w$ as one of your eigenvectors.
Given:
$$A = \begin{pmatrix} 3 &-4 &1 &0 \\ 4& 3 &0 &1 \\ 0 &0 &3 &-4 \\ 0 & 0 & 4 & 3 \end{pmatrix}$$
You have correctly found the eigenvalues:
$$ \lambda_1 = 3+4i, \lambda_2 = 3-4i$$
Each eigenvalue has an algebraic multiplicity of two.
For $\lambda_1 = 3 + 4i$, the row reduced echelon form (RREF) of $[A-\lambda_1 I]v_1 = 0$, is:
$$\begin{pmatrix} 1 &-i & 0 &0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}v_1 = 0$$
For this eigenvalue, you have correctly found one of the corresponding eigenvectors
$$v_1 = \begin{pmatrix} i \\ 1 \\ 0 \\ 0\\ \end{pmatrix}$$
However, we only have a single eigenvector (geometric multplicity is one) and need to find a second generalized eigenvector. One approach is to use $[A - \lambda_1 I]v_2 = v_1$, yielding the augmented matrix:
$$
\left[\begin{array}{rrrr|r}
4i &-4 &1 &0 & i \\ 4& 4i &0 &1 &1 \\ 0 &0 & 4i &-4 &0\\ 0 & 0 & 4 & 4i&0
\end{array}\right]
$$
The RREF of this matrix is:
$$
\left[\begin{array}{rrrr|r}
1 &-i & 0 & 0 & 0 \\ 0& 0 &1 &0 &i \\ 0 &0 & 0 &1 &1 \\ 0 & 0 & 0 & 0&0
\end{array}\right]
$$
From this RREF we have:
$$d = 1, c = i, a = i b$$
Choose $b = 0 \implies a = 0$, yielding a second generalized eigenvector:
$$v_2 = \begin{pmatrix} 0 \\ 0 \\ i \\ 1 \end{pmatrix}$$
It is worth noting that the eigenvalues come in complex conjugate pairs and so do the eigenvectors. In other words, you can just write the other two eigenvectors for $\lambda_2$ from the work above by taking the complex conjugate of each eigenvector.
Best Answer
If you just want to find all the eigenvectors corresponding to your chosen eigenvalue, then you can just find the basis for $ker(A-\lambda I)^n$ where $n$ is the multiplicity of the eigenvalue in the minimal polynomial. Because, all the lower power kernels are contained in the $n^{th}$ power kernel.
If you want to build Jordan chains for transformation matrix for Jordan canonical form, however, your chain vectors should obey the relation you described in your question. That is, for example if we have a sample case below:
Then your vectors in a chain should be related by
$f_i^k = (A-\lambda I)f_{i+1}^k$.
Also they should be linearly independent.
It is up to you how to choose those vectors, provided that they obey the above relations.
The chains and their length can be found by checking the dimensionality of discs (dimension difference of the consecutive Kernels). In other words, if there are two linearly independent vectors in the third disc, then you are going to have two chains with length 3. Then, after subtracting these two vectors, you can find only one linearly independent vector in the second disc and therefore you have one chain with length 2. And lastly, in the innermost circle you have one remaining linearly independent vector, so your last chain has length 1.
It can be easy to start from the outermost disc, because when you go from outer to inner disc, you just multiply by $(A-\lambda I)$. Also, note that when you choose independent vectors in some disc, and then multiply them by $(A-\lambda I)$ to get the corresponding vectors in the inner disc, those new vectors are going to be linearly independent. It also eases the procedure. But in that case finding $f_3^k$'s can be difficult, because besides being in the kernel of $(A-\lambda I)^3$ they should not be in the kernel of $(A-\lambda I)^2$.
When starting from innermost disc, it is easy to find normal eigenvectors. The vectors in the innermost disc are just normal eigenvectors (not generalized). However, then for building up your chain you need to find a vector in the outer disc so that $f_i^k = (A-\lambda I)f_{i+1}^k$ holds.