why is the PDF of the sum of two continuous random variables the convolution of the PDF's? I have taken a probability course last year but we didn't cover that notion. 
I do know the steps in the discrete case. finding the support of $X + Y$... calculating the the probability of each element of the support... dressing a table. however this process is very intuitive and self-explanatory.
the convolution is very useful when taking inverse of product of Laplace/Fourier transforms which is why it's hard for me to think of an analogy between taking the inverse of integral and computing the PDF of the sum of two random variables. 
I'd like to know the intuition behind : 
$$f_Z(z)=\int^{\infty}_{-\infty}f(x,z-x)dx$$
or just $$f_Z(z) = \int_{- \infty}^{\infty} f_X(x)f_Y(z-x)\;dx $$
when they are independent 
not a proof as I already found some on this site and elsewhere.
 A: I don't know what proofs you've seen, but the basic idea is that, if $X+Y=Z$, then $Y = Z-X$.  Then take all possible values for $X$, and add them up.
Let's pass into discrete probability, just to be simple.  Let $X$ and $Y$ be discrete random variables, and $Z = X+Y$.  Notice that
$$
    P(Z=z \mid X = x) = P(Y = z-x)
$$
So
$$
    P(Z=z) = \sum_{x} P(Z = z \mid X= x) P(X=x) = \sum_x P(Y=z-x) P(X=x)
$$
If you set $f_X(x) = P(X=x)$ and so on, this amounts to 
$$
    f_Z(z) = \sum_{x} f_X(x) f_Y(z-x)
$$
A: As a bit of graphical intuition:
Let's say you have two independent datasets, like this:
%pylab
x = randn(10000) + randint(0, 2, 10000) * 5
y = randn(10000) + randint(0, 2, 10000) * 5

fig, axes = plt.subplots(
    2, 2, sharex='col', sharey='row',
    gridspec_kw={'width_ratios': [1, 0.1], 'height_ratios': [0.1, 1]})

axes[1, 0].scatter(x, y, s=1)
axes[0, 0].hist(x, bins=200)
axes[1, 1].hist(y, bins=200, orientation='horizontal')

for z in linspace(-2, 8, num=6):
    # z = x + y  =>  y = z - x
    # which always has a slope of -1 for any fixed z.
    axes[1, 0].axline((z, z), slope=-1, c='orange')

Which look like this:

The orange lines are lines of equal z values.
A histogram of z is thus just the sum of the number of points along each of these lines (with width equal to the histogram bin width):
hist(x + y, bins=200)


In the continuous case, this is a convolution, as described in other answers. You're basically  asking:

For every possible fixed z value, what is the integral of all of the probabilities along that y = z - x line, relative to the entire joint PDF?

A: In the following I write the derivation using independent variables, but you can just substitute the product with the compound version. If you use the PDFs it becomes very apparent: Let $Z=X+Y$.
$$
\begin{aligned}
F_Z(t) = P(Z \leq t) = P(X+Y \leq t)
\end{aligned},
$$
$$
\begin{aligned}
F_Z(t) &= \iint\limits_{x+y\leq t}f(x,  y)dydx, \\
&= \int\limits_{-\infty}^\infty  \int\limits_{-\infty}^{t-x}f(x, y)dydx  \\
& \qquad \qquad\qquad \qquad  z := x + y \quad \text{(Substitution)}\\
&= \int\limits_{-\infty}^\infty  \int\limits_{-\infty}^{t}f(x, z-x)dzdx
\\
&= \int\limits_{-\infty}^{t} \int\limits_{-\infty}^\infty  f(x, z-x) dxdz \\
\end{aligned}
$$
Then
$$f_Z(t) = \frac{d}{dt} F_Z(t) = \int\limits_{-\infty}^\infty  f(x, t-x)dx$$
A: Think about what a density function represents in terms of probabilities in the discrete, independent case. 
$f_Z(z) = P(Z = z) = P(X+Y = z) = \sum_x P(X=x \,\,\&\,\, Y = z-x)=\sum_x f_X(x) f_Y(z-x)$. 
Your expression in the independent case is the continuous analog of this.
