$\newcommand{\D}{\mathcal{D}}$ In section 14.6 of Schwartz's "QFT and the Standard Model (7$\,^{\text{th}}$ printing)" [1], the author calculates the exact free fermionic partition function and, from that, derives the 2-point function (propagator). However, I am having trouble following the precise details of the latter derivation. I will state what's given in the text, then explicate my confusion.
The partition function for a free fermion (just by itself, i.e. not coupled to any gauge fields) is:
$$Z[\bar{\eta},\eta]=\int \D \bar{\psi}\D\psi\exp\left(i\int d^4x\left[\bar{\psi}(i\not\partial-m+i\epsilon)\psi+\bar{\eta}\psi+\bar{\psi}\eta\right]\right)\tag{14.101}$$
We can (formally) evaluate this explicitly.
$$Z[\bar{\eta},\eta]=N\exp\left(i\int d^4x \int d^4y \,\left[\bar{\eta}(x)(i\not\partial-m+i\epsilon)^{-1}\eta(y)\right]\right) \tag{14.102}$$
With this we can calculate the 2-point function.
$$\begin{align} \langle 0|T\left\{\psi(x)\bar{\psi}(y)\right\}|0\rangle&=(-i)^2\frac{1}{Z[0,0]}\left.\frac{\partial^2Z}{\partial \bar{\eta}(x)\partial \eta(y)}\right|_{\eta=\bar{\eta}=0}\\ &\overset{\color{red}{(?)}}{=}\frac{i}{i\not\partial-m+i\epsilon}\delta^{(4)}(x-y)\\ &=\int \frac{d^4p}{(2\pi)^4}\frac{i}{\not p-m+i\epsilon}e^{ip(x-y)} \tag{14.103} \end{align}$$
I have placed a red question mark over the step that confuses me.
[My Attempt:] The functional derivatives that I'll be taking are defined in the following way:
$$f(x)=\int d^4y \,G(x,y)g(y)\longrightarrow \frac{\partial f(x)}{\partial g(y)}=G(x,y) \tag{1}$$
which by the way immediately implies:
$$\frac{\partial f(x)}{\partial f(y)}=\delta^4(x-y)\tag{2}$$
So with this, I'll take the relevant functional derivatives of $Z[\bar{\eta},\eta]$.
$$\begin{align} \frac{\partial Z}{\partial \eta(y)}&=\frac{\partial}{\partial \eta(y)}\left(N\exp\left(i\int d^4x_1 \int d^4x_2 \,\left[\bar{\eta}(x_1)(i\not\partial-m+i\epsilon)^{-1}\eta(x_2)\right]\right) \right)\tag{3}\\ &=\left(i\int d^4 x_1 \bar{\eta}(x_1)(i\not\partial_y-m+i\epsilon)^{-1}\right)Z[\bar{\eta},\eta]\tag{4} \end{align}$$
Acting again gives us:
$$\begin{align} \frac{\partial^2Z}{\partial \bar{\eta}(x)\partial \eta(y)}&=(i(i\not\partial_y -m +i\epsilon)^{-1})Z[\bar{\eta},\eta]\\ &\,\,\,\,+\left(i\int d^4 x_1 \bar{\eta}(x_1)(i\not\partial_y-m+i\epsilon)^{-1}\right)\left(i\int d^4 x_2 (i\not\partial-m+i\epsilon)^{-1}\eta(x_2)\right)Z[\bar{\eta},\eta]\tag{5} \end{align}$$
Finally, setting $\eta=\bar{\eta}=0$ and dividing by $-Z[0,0]$ I get:
$$\begin{align}\langle 0|T\left\{\psi(x)\bar{\psi}(y)\right\}|0\rangle &=(-i)^2\frac{1}{Z[0]}\left.\frac{\partial^2Z}{\partial \bar{\eta}(x)\partial \eta(y)}\right|_{\eta=\bar{\eta}=0}\\ &=\boxed{\frac{-i}{i\not\partial_y-m+i\epsilon}}\tag{6} \end{align}$$
which is different from the author's (correct) answer. Where have I gone wrong? I am sure my confusion lies in taking functional derivatives.
[EDIT] I have read this related question, but I don't think anything there provides a direct answer to my question. That question concerns the precise meaning of inverting an operator, whereas mine is about taking functional derivatives properly and getting from them appropriate $\delta$-function factors.
