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 :)
Best Answer
Roughly speaking, you can interprete Wick theorem for quantum fields as the generalization of Wick theorem for Gaussian random variable, which also called Isserlis' theorem.
Omitting conceptual details, I can speak that $\phi(x_1)$, ... are Gassian distributed. According to the mentioned theorem, I can simply write the following: $$\langle\phi(x_1)\phi(x_2)\phi(y_1)^4\rangle = \text{sum of all possible two-point averages},$$ where the right hand side is given by (correct me if I miss something): $$\langle\phi(x_1)\phi(x_2)\rangle\langle\phi(y_1)\phi(y_1)\rangle\langle\phi(y_1)\phi(y_1)\rangle+\langle\phi(x_1)\phi(y_1)\rangle\langle\phi(y_1)\phi(y_1)\rangle\langle\phi(y_1)\phi(x_2)\rangle.$$ Next, you know that in representation picture you have integration over space, so you should add integration over $y_1$. From free theory analysis, you know that $D(x-y)=\langle\phi(x)\phi(y)\rangle$ is the simple quantity, which is called free propagator. Finally, you two-point function becomes $$\int d^4y_1\,D(x_1-x_2)D(y_1-y_1)D(y_1-y_1)+\int d^4y_1\,D(x_1-y_1)D(y_1-y_1)D(y_1-x_2).$$
The key point of using Feynman diagrams is do not to repeat this boring procedure of writing down all the possible Wick contractions. From the lagrangian you can immediately identify what is the bare vertex and what is the free theory propagator. It is enough to write down any possible correction for two-point correlator function