Today I took a linear algebra test which had the following problem:
(a) If $N\in {\rm Mat}(3,\Bbb C)$ is nilpotent, then the nilpotence index of $N$ is less or equal than $3$.
(b) Show that $X= {\rm Id}_3 + \frac{1}{n}N + \frac{1-n}{2n^2}N^2$ satisfies $X^n={\rm Id}_3 +N$ for each positive integer $n$.
(c) Conclude that every non-singular matrix $A\in {\rm Mat}(3,\Bbb C)$ has a $n$th root, that is, there is $B\in {\rm Mat}(3,\Bbb C)$ such that $B^n = A$.
I do not need to know how to solve it, since I managed to do it there. But I wondered if this generalizes for dimensions greater than $3$. For example, could one come up with some expression for $X$ in item (b) with some extra higher powers of $N$, so that repeating the straightforward reasoning using the Jordan canonical form in (c) works well?
Yep. Note that we can find a formal power series $f = \sum_{n=0}^{\infty} a_n X^n \in \mathbb{C}[[X]]$ such that $f^k = 1 + X$. Explicitly, the coefficients $a_n$ are given by $a_n = { \frac{1}{k} \choose n }$ which is consistent with the real/complex power series expansion of $(1+x)^{\frac{1}{k}}$ around $x = 0$ and the first three terms of $f$ are precisely the coefficients you have in part $(b)$.
Note that if $N \in M_l(\mathbb{C})$ is nilpotent, we can define a map $\operatorname{ev}_N \colon \mathbb{C}[[X]] \rightarrow M_l(\mathbb{C})$ by plugging in the matrix $N$ into a formal power series:
$$ \operatorname{ev}_N \left( \sum_{n=0}^{\infty} b_n X^n \right) = \sum_{n=0}^{\infty} b_n N^n $$
and this map is easily seen to be a ring homomorphism.
Given a nilpotent $N \in M_l(\mathbb{C})$ and $0 \neq \lambda \in \mathbb{C}$, define a matrix $B \in M_{l}(\mathbb{C})$ by
$$ B "=" (\lambda I_l + N)^{\frac{1}{k}} "=" \lambda^{\frac{1}{k}} (I_l + \frac{N}{\lambda})^{\frac{1}{k}} := \lambda^{\frac{1}{k}} \operatorname{ev}_{\frac{N}{\lambda}} (f) = \sum_{n=0}^{\infty} a_n \lambda^{\frac{1}{k}-n} N^n. $$
Since $\operatorname{ev}_N$ is a ring homomorphism, we have
$$ B^k = \lambda (\operatorname{ev}_{\frac{N}{\lambda}}(f))^k = \lambda \operatorname{ev}_{\frac{N}{\lambda}}(f^k) = \lambda \operatorname{ev}_{\frac{N}{\lambda}}(1 + X) = \lambda(I + \frac{N}{\lambda}) = \lambda I + N $$
and so $B$ is a $k$-th root of $\lambda I + N$. Then, like you wrote, by applying the Jordan decomposition we see that any invertible matrix has a $k$-th root.