Question: I would like to know why an polynomial fit is divergent on the boundaries of the fit interval (as shown on the doc site of the matlab polyfit function).
And why a polynomial fit on a data set is a very bad extrapolation technique for data just outside of the fit interval.

Method used to compute coefficients: least square method.

It seems that splines are better techniques for that but I try to understand why polynomial fits have these drawbacks.

Answer: 
Yes! Pick your favorite variable, say $a$, and treat the expression as a polynomial in $a$ (with coefficients that are themselves polynomials in $b,c,d$). It's easy to check that this polynomial has degree at most $2$; but it's also easy to check that it vanishes at $b$, $c$, and $d$ and thus has at least three roots. Therefore it must be the zero polynomial.
Making two variables equal makes two products vanish and makes the remaining two cancel each other. E.g. with $a=b$, the equation turns to the identity

$$(c-d)(b-d)(b-c)-(c-d)(b-d)(b-c)\equiv 0.$$ Hence the given expression is of the form

$$P(a,b,c,d)(a-b)(a-c)(a-d)(b-c)(b-d)(c-d).$$

Because of the degrees, $P(a,b,c,d)\equiv 0$.