From L. Evans' book on PDE's:
Definition
The function
$$\Phi(x) = \begin{cases} -\frac{1}{2\pi} \log|x| & (n=2) \\ \frac{1}{n(n-2)\alpha(n)} \frac{1}{\lvert x \rvert^{n-2}} & (n\geq3)\end{cases}$$
defined for $x\in\mathbb{R}^n, \ x \neq 0$, is the fundamental solution of Laplace's equation.
By construction, the function $x \mapsto \Phi(x)$ is harmonic for $x\neq0$. If we shift the origin to a new point $y$, the PDE $\Delta u=0$ is unchanged, and so $x\mapsto\Phi(x-y)$ is also harmonic as a function of $x$, $x\neq y$. Let us now take a $f:\mathbb{R}^n\rightarrow\mathbb{R}$ and note that the mapping $x\mapsto\Phi(x-y)f(y)$, $x\neq y$ is harmonic for each $y\in\mathbb{R}^n$, and thus so is the sum of finitely many such expressions built for different points $y$.
Now I know that a function is harmonic if it is $\mathbf{C}^2$ and it satisfies $\Delta u=0$, and that a set of harmonic functions can be seen as a vector space, i.e. linear combinations of harmonic functions are again harmonic functions, but I don't know how to show the function $x\mapsto\Phi(x)$ is harmonic, nor how to show that the mapping $x\mapsto\Phi(x-y)f(y)$, $x\neq y$ is harmonic for each $y\in\mathbb{R}^n$, when $x\mapsto\Phi(x)$ is harmonic.