Actually calculating something using Wick's Theorem I am still struggling to get my head around QFT and whilst I think I understand the method of generating functionals to compute correlation functions (as in my question here), my course notes often use Wick's theorem for calculations, and I'd really like to understand how that works.
So if I return to the standard example of scalar Yukawa theory with interaction Lagrangian:
$$\mathcal{L}_{\text{int}}=-g\psi^{\dagger}\psi\phi$$
I would like to know how to compute quantities like:
$$\langle \Omega|T[\psi^{\dagger}(x)\psi(y)\phi(z)]|\Omega\rangle$$
Where $T[...]$ represents the time-ordering operator. My understanding is that this can be expressed as:
$$\langle \Omega|T[\psi^{\dagger}(x)\psi(y)\phi(z)]|\Omega\rangle = \lim_{t\to\infty}\frac{\langle0|T[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)\exp\left(-i\int_{-t}^{t}\mathcal{H}_{\text{int}}\:\mathrm{d}^{4}s\right)]|0\rangle}{\langle0|T[\exp\left(-i\int_{-t}^{t}\mathcal{H}_{\text{int}}\:\mathrm{d}^{4}s\right)]|0\rangle}$$
We can expand the exponential to (arbitrarily) first-order in $g$:
$$\langle0|T[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)\exp\left(-i\int_{-t}^{t}\mathcal{H}_{\text{int}}\:\mathrm{d}^{4}s\right)]|0\rangle = \langle 0|T\left[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)(1-ig\int_{-t}^{t}\psi_{I}^{\dagger}(s)\psi_{I}(s)\phi_{I}(s)\:\mathrm{d}^{4}s)\right]|0\rangle$$
By Wick's theorem, the zeroth order in $g$ has only 3 operators and so will inevitably yield a contribution of 0. The second gives:
$$\langle0|\int_{-t}^{t}T[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)\psi_{I}^{\dagger}(s)\psi_{I}(s)\phi_{}(s)]\:\mathrm{d}^{4}s|0\rangle$$
However, I'm not sure what the best way to proceed is. I can use Wick's theorem to expand the time-ordered string of operators in terms of products of operators, but then where do I go?
 A: In your example there's two contractions giving two terms. Its amplitude is
$$ \left< \psi^{\dagger}_a(x) \psi^b(y) \phi(z) \right> = -i \lambda \int d^4 s \, \Delta_M(z - s) \Delta_m(y-s)^b_{\;c} \Delta_m(x-s)^c_{\;a} $$
$$ - i \lambda \int d^4 s \Delta_M(x - y) ^b_{\;a} \Delta_M(0)^c_{\;c} \Delta_m(s-x) + \mathcal{O}\left(\lambda^2\right). $$
How can I see that this is true? The algorithm is pretty simple, actually:


*

*You have a term corresponding to each possible contraction. A contraction is a diagram where pairs of $(\psi^{\dagger}, \psi)$ and $(\phi, \phi)$ are replaced by the corresponding propagators.

*You have to exclude the diagrams containing bubble subgraphs. A bubble graph is a graph with no external legs. An external leg is a contraction which has one of the fields appearing in the correlation bracket ($\psi^{\dagger}(x), \psi(y), \phi(z)$). This is because we want to account for the normalization $\mathcal{N}/\mathcal{N}_0$ as I mentioned in the answer to your previous question. Proofs of this can be found in any QFT textbook, e.g. Peskin-Schreder.

*For each internal (interaction) vertex we have a factor of $$-i \lambda \int d^4 s. $$

*Each term is a product of integrals over spacetime positions of internal (interaction) vertices and propagators.
In your example there's two $\phi$ fields, two $\psi$ fields and two $\psi^{\dagger}$ fields. Therefore we must contract $\phi$ with $\phi$, but we have a choice of which $\psi^{\dagger}$ gets contracted with specific $\psi$.
These correspond to the two diagrams below:


Both don't contain bubble subgraphs, however the second contains the tadpole divergence $\Delta_m(0)$. These divergences arise in QFT often. They have to be renormalized by requiring that $\left< \phi \right> = 0$, which is eqiuvalent to throwing away the tadpole contribution. (Actually, we don't just throw away mathematical expressions; we absorb them in the redefinition of fields).
The relevant part of your vertex amplitude is thus
$$ \left< \psi^{\dagger}_a(x) \psi^b(y) \phi(z) \right> = -i \lambda \int d^4 s \, \Delta_M(z - s) \Delta_m(y-s)^b_{\;c} \Delta_m(x-s)^c_{\;a}. $$
Let me know if you have further questions.
P.S. oh and I supposed that $\psi$ and $\psi^{\dagger}$ are field multiplets, thus the internal indices $a, b, c$ labeling the components of the multiplet. If it's just a complex number then just throw these indices away :)
A: I myself had the same question once. There is not much philosophy. Just perform the contractions. 
Wick's theorem states that only fully contracted terms contribute to correlation functions. So, you have to start contracting things until there is nothing else.
This is a boring process, and it is easier with Feynman diagrams. Allow me to do another problem, which is faster to type.
Warning : I am ignoring factors of i, $\frac{1}{2}$, etc. Just follow the logic.
Consider the theory with $H_{\text{int}} = \lambda \int d^{3}x \, \phi^{4}(x)$ and let us calculate the 2-point function to first order in $\lambda$: 
$$
\langle 0 | T \phi(x) \phi(y) \int dt \, H_{\text{int}} | 0 \rangle
$$
Let us write it as follows
$$
\langle 0 | T \phi(x) \phi(y) \int d^{4} z \,\, \phi(z)\,\phi(z)\,\phi(z)\,\phi(z) | 0 \rangle.
$$
Wick's theorem states that only full contractions survive. This means that there must not remain a single $\phi(x)$ floating around.
In this example, there are two independent Wick contractions:


*

*Contract $\phi(x)$ with $\phi(y)$ and $\phi(z)$'s between themselves;

*or contract $\phi(x)$ with $\phi(z)$, $\phi(y)$ with $\phi(z)$ and $\phi(z)$ with $\phi(z)$.
The first contraction can be done in three ways. Since the $x$ and $y$ contraction is fixed, you just have to choose a partner for $\phi(z)$, and there are $3$ of them. The second contraction can be done in 12 ways. There are 4 choices to contract $\phi(x)$ with $\phi(z)$ and then $3$ choices remain when you contract $\phi(y)$ with $\phi(z)$.
The result is then
$$
\langle 0 | T \phi(x) \phi(y) \int d^{4} z \,\, \phi(z)\,\phi(z)\,\phi(z)\,\phi(z) | 0 \rangle = 3 \lambda \Delta(x-y) \int d^{4} z \, \Delta (z-z) \, \Delta(z-z) \\
+ 12 \lambda \int d^{4} z \, \Delta(x-z) \, \Delta(y-z) \Delta(z-z).
$$
You could now work at order $\lambda^{2}$: 
$$
\langle 0 | T \phi(x) \phi(y) \int d^{4} z \,\, \phi(z)\,\phi(z)\,\phi(z)\,\phi(z)  \int d^{4} u \,\, \phi(u)\,\phi(u)\,\phi(u)\,\phi(u)| 0 \rangle.
$$
But it is quickly recognizable how this is going to be a pain in the butt. A way around is to use Feynman diagrams.
In a diagram, a contraction is represented by a line and the $H_{\text{int}}$ becomes what is known as a vertex. You can see how it is a vertex, because it is precisely in the interacting Hamiltonian where the non-linearities are. So all contractions with the terms inside the Hamiltonian will give raise to vertices of the diagram.
I would draw the diagram for the example that I just did, but I do not know how to do this in here. Sorry.
