If I interpret your initialization statement correctly, you have
Y(a,b,0,0) = 1
for every non-negative integers a and b. You can take advantage of this in evaluating Y for arbitrarily large b.
As to the parameter a, it doesn't change within the Y's in the recursion formula, so you can consider this as a separate recursion problem for each different a.
I would recommend the following ordering in calling on the recursion formula:
Do this for each a:
for i = 1:N
for j = 1:i
b = i-j;
for k = 0:j % <-- Corrected. The k = 1:j was wrong
c = j-k;
d = k;
Y(a,b,c,d) = 1/(c+d)*(.....)
This way, each of the Y values needed to compute Y(a,b,c,d) will already have been computed provided you make use of your initialization requirements. On exit from the for-loops Y will have been computed for each non-negative integer value of b, c, and d such that b+c+d<=N.
