I don't know where else to find this result, but it is not hard to prove. Indeed, without loss of generality $\bar x=0$. So, we have $f(0)=0$, $M^n=0$ for $M:=f'(0)$, and $\|f'(x)-M\|\le2L\|x\|$ for some real $L\ge0$, some norms on $\mathbb R^{n\times n}$ and $\mathbb R^n$ (both denoted here by $\|\cdot\|$), some real $h>0$, and all $x\in B_h:=\{x\in\mathbb R^n\colon\|x\|\le h\}$. So, for $x\in B_h$ and
$$g(x):=f(x)-Mx\tag{1}$$
one has
$$\|g(x)\|=\|f(x)-f(0)-Mx\|=\Big\|\int_0^1(f'(tx)-M)x\,dt\Big\| $$
$$
\le\int_0^1\|f'(tx)-M\|\,\|x\|\,dt\le\int_0^1 2Lt\|x\|\,\|x\|\,dt=L\|x\|^2,
$$
whence $\|f(x)\|\le\|M\|\|x\|+L\|x\|^2\le(\|M\|+Lh)\|x\|$.
In fact, letting $\|M\|$ be the (say) Euclidean operator norm of the Jordan form of the nilpotent matrix $M$, we may assume that $\|M\|\le1$, so that $\|f(x)\|\le(1+Lh)\|x\|$ for $x\in B_h$. Hence, $x_j\in B_h$ for all $x\in B_{h_1}$ and $j=0,\dots,n-1$, where $h_1:=h/(1+Lh)^{n-1}$, and, for each $x$, the sequence $(x_j)$ is defined recursively by the formulas $x_0:=x$ and $x_{j+1}:=f(x_j)$.
By $(1)$, $f(x)=Mx+g(x)$. So,
in view of the condition $M^n=0$, by induction, for any natural $N\ge n-1$ one obtains
$$x_N=M^Nx+M^{N-1}g(x)+\dots+Mg(x_{N-2})+g(x_{N-1})$$
$$
=M^{n-1}g(x_{N-n})+\dots+Mg(x_{N-2})+g(x_{N-1})
$$
for $x\in B_{h_2}$, where $h_2:=\min[h_1,\frac1{nL}]$;
here we use the estimate
$$\|x_N\|\le\|M^{n-1}g(x_{N-n})+\dots+Mg(x_{N-2})+g(x_{N-1})\|$$
$$
\le\|g(x_{N-n})\|+\dots+\|g(x_{N-2})\|+\|g(x_{N-1})\|
$$
$$
\le L\|x_{N-n}\|^2+\dots+L\|x_{N-1}\|^2\le nLh_2^2\le h_2
$$
for $x\in B_{h_2}$. In particular, it follows that
$$\|x_N\|
\le L(\|x_{N-n}\|^2+\dots+\|x_{N-1}\|^2)
$$
for $x\in B_{h_2}$, which
yields the desired quadratic convergence.
