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.
Write it as$$\begin{align}&\int\_a^b(c-d)[2x-c-d]\mathrm{d}x+\int\_c^d(a-b)[2y-a-b]\mathrm{d}y\\&=\int\_a^b\mathrm{d}x\int\_c^d\mathrm{d}y\{2x+2y-a-b-c-d\}.\end{align}$$The symmetry $x\to(a+b) /2-x,\,y\to(c+d)/2-y$ make the rest obvious.