The analytical result for free massless fermion propagator

Method One:

\begin{eqnarray*} & & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x}\\ & = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\int dk_{0}e^{-ik_{0}x_{0}}\frac{1}{[k_{0}+(|\mathbf{k}|-i\epsilon)][k_{0}-(|\mathbf{k}|-i\epsilon)]}\\ & = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\bigg(\theta(x_{0})(-2\pi i)\frac{1}{2|\mathbf{k}|}e^{-i|\mathbf{k}|x_{0}}+\theta(-x_{0})(2\pi i)\frac{1}{-2|\mathbf{k}|}e^{i|\mathbf{k}|x_{0}}\bigg)\\ & = & \frac{\pi}{(2\pi)^{4}}\int d^{3}k\frac{1}{|\mathbf{k}|}e^{i\mathbf{k}\cdot\mathbf{x}}e^{-i|\mathbf{k}||x_{0}|}\\ & = & \frac{\pi}{(2\pi)^{4}}2\pi\int_{0}^{\infty}d|\mathbf{k}||\mathbf{k}|^{2}\frac{1}{|\mathbf{k}|}e^{-i|\mathbf{k}||x_{0}|}\int_{-1}^{1}dye^{i|\mathbf{k}||\mathbf{x}|y}\\ & = & \frac{\pi}{(2\pi)^{4}}2\pi\int_{0}^{\infty}d|\mathbf{k}||\mathbf{k}|e^{-i|\mathbf{k}||x_{0}|}\frac{1}{i|\mathbf{k}||\mathbf{x}|}(e^{i|\mathbf{k}||\mathbf{x}|}-e^{-i|\mathbf{k}||\mathbf{x}|})\\ & = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\int_{0}^{\infty}d|\mathbf{k}|(e^{-i|\mathbf{k}|(|x_{0}|-|\mathbf{x}|)}-e^{-i|\mathbf{k}|(|x_{0}|+|\mathbf{x}|)})\\ & = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\bigg[\bigg(\pi\delta(|x_{0}|-|\mathbf{x}|)-\mathscr{P}\frac{i}{|x_{0}|-|\mathbf{x}|}\bigg)-\bigg(\pi\delta(|x_{0}|+|\mathbf{x}|)-\mathscr{P}\frac{i}{|x_{0}|+|\mathbf{x}|}]\bigg)\bigg]\\ & = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\bigg[\pi\delta(|x_{0}|-|\mathbf{x}|)+\mathscr{P}\bigg(\frac{i}{|x_{0}|+|\mathbf{x}|}-\frac{i}{|x_{0}|-|\mathbf{x}|}\bigg)\bigg]\\ & = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}\bigg(2\pi|\mathbf{x}|\delta(x^{2})-2i|\mathbf{x}|\mathscr{P}\frac{1}{x^{2}}\bigg)\\ & = & \frac{\pi}{(2\pi)^{4}}2\pi\frac{1}{i|\mathbf{x}|}(-2i|\mathbf{x}|)\bigg(\mathscr{P}\frac{1}{x^{2}}+i\pi\delta(x^{2})\bigg)\\ & = & -\frac{1}{4\pi^{2}}\frac{1}{x^{2}-i\epsilon} \end{eqnarray*}

i.e. \begin{eqnarray*} \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x} & = & -\frac{1}{(2\pi)^{2}}\frac{1}{x^{2}-i\epsilon} \end{eqnarray*}

where we have used \begin{eqnarray*} \frac{1}{x+i\epsilon} & = & \mathscr{P}\frac{1}{x}-i\pi\delta(x) \end{eqnarray*}

\begin{eqnarray*} \int_{0}^{\infty}e^{ikx}dk & = & \lim_{\epsilon\to0^{+}}\int_{0}^{\infty}e^{ik(x+i\epsilon)}dk=\lim_{\epsilon\to0^{+}}\frac{i}{x+i\epsilon}=\pi\delta(x)+\mathscr{P}\frac{i}{x} \end{eqnarray*} \begin{eqnarray*} \delta(x^{2}) & = & \delta(x_{0}^{2}-\mathbf{x}^{2})=\frac{1}{2|\mathbf{x}|}[\delta(x_{0}-|\mathbf{x}|)+\delta(x_{0}+|\mathbf{x}|)]\\ & = & \frac{1}{2|\mathbf{x}|}[\theta(x_{0})\delta(x_{0}-|\mathbf{x}|)+\theta(-x_{0})\delta(x_{0}+|\mathbf{x}|)]\\ & = & \frac{1}{2|\mathbf{x}|}\delta(|x_{0}|-|\mathbf{x}|) \end{eqnarray*}

Method Two:

We can also see that \begin{eqnarray*} & & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+m^{2}+i\epsilon}e^{-ik\cdot x}\\ & = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\int dk_{0}e^{-ik_{0}x_{0}}\frac{1}{[k_{0}+(E_{k}-i\epsilon)][k_{0}-(E_{k}-i\epsilon)]}\\ & = & \frac{i}{(2\pi)^{4}}\int d^{3}ke^{i\mathbf{k}\cdot\mathbf{x}}\bigg(\theta(x_{0})(-2\pi i)\frac{1}{2E_{k}}e^{-iE_{k}x_{0}}+\theta(-x_{0})(2\pi i)\frac{1}{-2E_{k}}e^{iE_{k}x_{0}}\bigg)\\ & = & \frac{i(-2\pi i)}{2(2\pi)^{4}}\int d^{3}k\frac{1}{E_{k}}e^{i\mathbf{k}\cdot\mathbf{x}}e^{-iE_{k}|x_{0}|}\\ & = & \frac{i(-2\pi i)(2\pi)}{2(2\pi)^{4}}\int_{0}^{\infty}dkk^{2}\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\int_{-1}^{1}dye^{ik|\mathbf{x}|y}\\ & = & \frac{i(-2\pi i)(2\pi)}{2(2\pi)^{4}}\int_{0}^{\infty}dkk^{2}\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{1}{ik|\mathbf{x}|}(e^{ik|\mathbf{x}|}-e^{-ik|\mathbf{x}|})\\ & = & \frac{i(-2\pi i)(2\pi)2}{2(2\pi)^{4}}\int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}\\ & = & \frac{1}{(2\pi)^{2}}\int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}\\ & \equiv & \frac{1}{(2\pi)^{2}}\bigg(\theta(x^{2})\times\mathrm{I}+\theta(-x^{2})\times\mathrm{II}\bigg) \end{eqnarray*}

If $x^{2}>0$, we can choose a frame with $x^{\mu}=(x_{0},\mathbf{0})$ and $x^{2}=x_{0}^{2}$. We can see \begin{eqnarray*} \mathrm{I} & = & \int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}=\int_{0}^{\infty}dkk^{2}\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\\ & = & \int_{m}^{\infty}dE_{k}\sqrt{E_{k}^{2}-m^{2}}e^{-iE_{k}|x_{0}|}\\ & = & m^{2}\int_{1}^{\infty}dt\sqrt{t^{2}-1}e^{-im|x_{0}|t},\ \bigg[\text{note: }a\equiv m|x_{0}|=m\sqrt{x^{2}}\bigg]\\ & = & m^{2}\int_{1}^{\infty}dt(t^{2}-1)\frac{e^{-iat}}{\sqrt{t^{2}-1}}\\ & = & -m^{2}(\frac{\partial^{2}}{\partial a^{2}}+1)\int_{1}^{\infty}dt\frac{e^{-iat}}{\sqrt{t^{2}-1}}\\ & = & \frac{i\pi m^{2}}{2}[H_{0}^{(2)\prime\prime}(a)+H_{0}^{(2)}(a)] \end{eqnarray*}

where we have used \begin{eqnarray*} \int_{1}^{\infty}dt\frac{e^{-iat}}{\sqrt{t^{2}-1}} & = & -\frac{i\pi}{2}J_{0}(a)-\frac{\pi}{2}N_{0}(a)=-\frac{i\pi}{2}H_{0}^{(2)}(a),\ (a>0) \end{eqnarray*}

From

\begin{eqnarray*} Z_{\nu}^{\prime} & = & Z_{\nu-1}-\frac{\nu}{x}Z_{\nu}\\ Z_{\nu}^{\prime} & = & -Z_{\nu+1}+\frac{\nu}{x}Z_{\nu} \end{eqnarray*}

we can see \begin{eqnarray*} Z_{0}^{\prime} & = & -Z_{1}\\ Z_{0}^{\prime\prime} & = & -Z_{1}^{\prime}=-(Z_{0}-\frac{1}{x}Z_{1}) \end{eqnarray*}

i.e. \begin{eqnarray*} Z_{0}^{\prime\prime}+Z_{0} & = & \frac{1}{x}Z_{1} \end{eqnarray*}

So we have \begin{eqnarray*} \mathrm{I} & = & \frac{i\pi m^{2}}{2}[H_{0}^{(2)\prime\prime}(a)+H_{0}^{(2)}(a)]=\frac{i\pi m^{2}}{2a}H_{1}^{(2)}(a)=\frac{i\pi m}{2\sqrt{x^{2}}}H_{1}^{(2)}(m\sqrt{x^{2}}) \end{eqnarray*}

If $x^{2}<0$, we can choose a frame with $x^{\mu}=(0,\mathbf{x})$ and $x^{2}=-|\mathbf{x}|^{2}$. We can see \begin{eqnarray*} \mathrm{II} & = & \int_{0}^{\infty}dkk\frac{1}{E_{k}}e^{-iE_{k}|x_{0}|}\frac{\mathrm{sin}(k|\mathbf{x}|)}{|\mathbf{x}|}=\frac{1}{|\mathbf{x}|}\int_{0}^{\infty}dkk\frac{\mathrm{sin}(k|\mathbf{x}|)}{\sqrt{m^{2}+k^{2}}}\\ & = & \frac{m}{|\mathbf{x}|}\int_{0}^{\infty}dt\frac{t\mathrm{sin}(m|\mathbf{x}|t)}{\sqrt{1+t^{2}}},\ \bigg[\text{note: }b\equiv m|\mathbf{x}|=m\sqrt{-x^{2}}\bigg]\\ & = & \frac{m^{2}}{b}\int_{0}^{\infty}dt\frac{t\mathrm{sin}(bt)}{\sqrt{1+t^{2}}}\\ & = & -\frac{m^{2}}{b}\frac{\partial}{\partial b}\int_{0}^{\infty}dt\frac{\mathrm{cos}(bt)}{\sqrt{1+t^{2}}}\\ & = & -\frac{m^{2}}{b}K_{0}^{\prime}(b) \end{eqnarray*}

where we have used \begin{eqnarray*} \int_{0}^{\infty}dt\frac{\mathrm{cos}(bt)}{\sqrt{1+t^{2}}} & = & K_{0}(b),\ (x>0) \end{eqnarray*}

From

\begin{eqnarray*} Z_{\nu}^{\prime} & = & Z_{\nu-1}-\frac{\nu}{x}Z_{\nu}\\ Z_{\nu}^{\prime} & = & -Z_{\nu+1}+\frac{\nu}{x}Z_{\nu} \end{eqnarray*}

we can see \begin{eqnarray*} Z_{0}^{\prime} & = & -Z_{1} \end{eqnarray*}

So we get \begin{eqnarray*} \mathrm{II} & = & -\frac{m^{2}}{b}K_{0}^{\prime}(b)=\frac{m^{2}}{b}K_{1}(b)=\frac{m}{\sqrt{-x^{2}}}K_{1}(m\sqrt{-x^{2}}) \end{eqnarray*}

Finally we have \begin{eqnarray*} \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+m^{2}+i\epsilon}e^{-ik\cdot x} & = & \frac{1}{(2\pi)^{2}}\bigg(\theta(x^{2})\times\frac{i\pi m}{2\sqrt{x^{2}}}H_{1}^{(2)}(m\sqrt{x^{2}})+\theta(-x^{2})\times\frac{m}{\sqrt{-x^{2}}}K_{1}(m\sqrt{-x^{2}})\bigg)\\ & = & \frac{m}{(2\pi)^{2}\sqrt{|x^{2}|}}\bigg(\theta(x^{2})\frac{i\pi}{2}H_{1}^{(2)}(m\sqrt{x^{2}})+\theta(-x^{2})K_{1}(m\sqrt{-x^{2}})\bigg) \end{eqnarray*}

then we have \begin{eqnarray*} & & \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x}\\ & = & \lim_{m\to0}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+m^{2}+i\epsilon}e^{-ik\cdot x}\\ & = & \lim_{m\to0}\frac{m}{(2\pi)^{2}\sqrt{|x^{2}|}}\bigg(\theta(x^{2})\frac{i\pi}{2}\frac{2i}{\pi m\sqrt{x^{2}}}+\theta(-x^{2})\frac{1}{m\sqrt{-x^{2}}}\bigg)\\ & = & \frac{1}{(2\pi)^{2}}\bigg(-\theta(x^{2})\frac{1}{x^{2}}+\theta(-x^{2})\frac{1}{-x^{2}}\bigg)\\ & = & -\frac{1}{(2\pi)^{2}}\frac{1}{x^{2}} \end{eqnarray*}

where we have used \begin{eqnarray*} H_{1}^{(2)}(x) & = & J_{1}(x)-iN_{1}(x)\rightarrow\frac{2i}{\pi x}\ \ \text{as}\ \ x\rightarrow0\\ K_{1}(x) & \rightarrow & \frac{1}{x}\ \ \text{as}\ \ x\rightarrow0 \end{eqnarray*}

So we get \begin{eqnarray*} \int\frac{d^{4}k}{(2\pi)^{4}}\frac{i}{k^{2}+i\epsilon}e^{-ik\cdot x} & = & -\frac{1}{(2\pi)^{2}}\frac{1}{x^{2}} \end{eqnarray*}

with $\frac{1}{x^{2}-i\epsilon}$ replaced by $\frac{1}{x^{2}}$.

Yes it is correct. The derivation in P&S is straightforward but I will expand on it a bit. The key observation is that \begin{equation} \int\frac{d^4k}{(2\pi)^4}e^{-ik\cdot(y-z)}\frac{i\gamma^{\mu}k_{\mu}}{k^2+i\epsilon} =-\gamma^{\mu}\partial_{\mu}\int\frac{d^4k}{(2\pi)^4}\frac{1}{k^2+i\epsilon}e^{-ik\cdot(y-z)}, \end{equation} where the integral on the right hand side is the Feynman propagator for a massless scalar. Performing the $k$ integrals to get to position space yields \begin{equation} \int\frac{d^4k}{(2\pi)^4}\frac{1}{k^2+i\epsilon}e^{-ik\cdot(y-z)}=\frac{i}{4\pi^2}\frac{1}{(y-z)^2-i\epsilon}. \end{equation} If you aren't sure about this last step, it is easier to consider the massive case and then take the limit as $m\rightarrow 0$ at the end. Schwinger parameters are also helpful for proving this identity.

Once you have transformed to position space you simply act with $-\gamma^{\mu}\partial_{\mu}$ to arrive at the final expression.

As alluded to in the other answer here, the integral can basically be evaluated in several ways. One of them is the one that OP has himself followed. (PS - OP, congratulations on completing that feat!)

Let me present here a way of computing this using Schwinger parameterization. We will use $$ \frac{1}{a} = \int_0^\infty d\tau e^{- \tau a} ~, a > 0~. $$ We want to compute $$ I = \int \frac{d^4k}{(2\pi)^4} \frac{i e^{- i k \cdot x}}{k^2+ i \epsilon} = \int \frac{dk^0}{2\pi} \frac{d^3k}{(2\pi)^3} \frac{i e^{- i k^0 t + i \vec{k} \cdot \vec{x} }}{(k^0)^2 - \vec{k}^2+ i \epsilon} $$ Do a Wick rotation $k^0 \to i k^0_E$, $t \to - i t_E$. Then, $$ I = \int \frac{dk_E^0}{2\pi} \frac{d^3k}{(2\pi)^3} \frac{ e^{ k_E^0 t + i \vec{k} \cdot \vec{x} }}{(k_E^0)^2 + \vec{k}^2+ i \epsilon} = \int \frac{d^4k}{(2\pi)^3} \frac{ e^{ - i k_E^0 t_E + i \vec{k} \cdot \vec{x} }}{ k^2 } $$ where in the last equation, we now have a Euclidean $k^2$ that is always positive over the range of integration. Now, we may use the Schwinger parameterization so that $$ I = \int_0^\infty d\tau \int \frac{d^4k}{(2\pi)^4} e^{- k^2 \tau - i k_E^0 t_E + i \vec{k} \cdot \vec{x}} $$ Now, we can do the integral of $k$ quite easily since $k^2 = \sum_i k_i^2$. This gives $$ I = \int_0^\infty d\tau \frac{1}{(4\pi)^2} \frac{e^{ - \frac{1}{4\tau} ( t_E^2 + \vec{x}^2 ) } }{\tau^2} $$ Now to perform the integral over $\tau$, define new integration variable $y = \frac{1}{4\tau} $. Then $$ I = \int_0^\infty dy \frac{1}{(2\pi)^2} e^{ - y ( t_E^2 + \vec{x}^2 ) } $$ This last integral is again the Schwinger parameter one. It converges and is nice so we compute it and find $$ I = \frac{1}{4\pi^2 ( t_E^2 + \vec{x}^2 ) } = - \frac{1}{4\pi^2 ( t^2 - \vec{x}^2 ) } = - \frac{1}{4\pi^2 x^2 } ~. $$ where in the last step, we have performed the inverse Wick rotation to go back to Lorentzian time.