Double sum with different lengths of the operators

Hi guys. I have a simple case of a double sum where the first sum is running from 1..k, the second from 1..l, where k,l might have different sizes. The sum refers to factorial of (i+j):

sum_(i=1)^k sum_(j=1)^l (i+j)!

Ok, so here I define the sum:

f=@(k,l) sum( factorial((1:k)+(1:l)) )

x=1:5

y=1:4 %say, I want to evaluate the sum where the first index runs utp 5 and the second up to 4 Then:

arrayfun(f,x,y)

What am I doing wrong? Maybe I have to define differently the matrix (x,y) which has to be assigned to the sum?

Many thanks!