Before interpreting them in more advanced language like "string diagrams" or "monoidal closed categories" it might be good to stress that Feynman diagrams are very elementary combinatorial objects which encode **contractions of tensors**.
By tensor I mean an array of numbers like $A=(A_{a,b,c})$ with indices $a,b,c$ running over some finite sets which need to be specified.
If you have such objects say $B_{abcd}$ and $C_{ab}$ you can construct new ones like
$$
T_{a,b,c,d}:=\sum_{e,f,\ldots,l} C_{ae}C_{bg} B_{eghf} C_{fi}C_{hk}B_{iklj}C_{jc}C_{ld}
$$
Obviously one can produce tons of similar examples of increasing complexity and it is desirable to have a way of encoding precisely such constructions. A natural way of doing that is basically to use pictures or graphs which is what Feynman diagrams are.

Linear algebra "done wrong", i.e., *matrix algebra* is the particular case where tensors have one (vectors) or two indices (matrices) only. Although, the $n$-dimensional determinant introduces a higher (Levi-Civita) tensor
$\epsilon_{i_1\ldots i_n}$ given by the sign of the permutation $i_1\ldots i_n$ (and zero if indices are repeated).

A significant part of functional analysis is about studying what happens when discrete summation
indices $a,b,\ldots$ become continuous variables that are integrated over. Then, matrices $C_{a,b}$ become kernels $C(x,y)$ which one can make sense of, e.g, as distributions, by invoking the Schwartz Kernel Theorem.

Feynman diagrams, as they are used in quantum field theory, typically correspond
to this infinite-dimensional generalization. For instance, the diagram for the expression $T_{abcd}$ above becomes the main second order contribution to the four-point function of the $\phi^4$ model in $d$ dimensions if one decides that the tensor $C_{ab}$
now becomes the kernel $C(x,y)$ of the operator $-\Delta+m^2 I$ in $\mathbb{R}^d$
and one lets the tensor $B_{abcd}$ become the kernel $B(x,y,z,u)$ of the distribution in $S'(\mathbb{R}^{4d})$
given by the action
$$
f\longmapsto\ \int_{\mathbb{R}^d} f(x,x,x,x)\ d^dx
$$
on test functions $f\in S(\mathbb{R}^{4d})$.

As for why this should have to do with the Laplace/stationary phase method, the reason is because Gaussian integration is an "algebraic" operation. Namely, it can be expressed as a differential operator (albeit of infinite order).
For example if $\mu$ is the centered Gaussian measure on $\mathbb{R}^n$ with covariance $C_{a,b}$, then for any polynomial
$P\in \mathbb{R}[x_1,\ldots,x_n]$,
$$
\int_{\mathbb{R}^n} P(x)\ d\mu(x)=\left.\exp\left(\frac{1}{2}\sum_{a,b=1}^n C_{a,b} \frac{\partial^2}{\partial x_a\partial x_b}\right)\ P(x)\ \right|_{x=0}\ .
$$

Note that Haar integration on $SU(n)$ can also be expressed as an infinite order differential operator (see my two answers to How to constructively/combinatorially prove Schur-Weyl duality? ).
So Feynman diagrams also appear in invariant theory/representation theory (see my answer to Who invented diagrammatic algebra? for some examples and pictures by Kempe in the case of the invariants of the binary quintic that are given explicitly in nondiagrammatic fashion in my answer to
Explicit formulas for invariants of binary quintic forms ).

4more comments