function s = binomial(x,y,n)
% Calculates explicitly the binomial expansion of (x+y)^n
s = 0;
for k=0:n
    an = x^(n-k) * y^k * factorial(n) / ( factorial(k) * factorial(n-k) );
    s = s + an;
end
