I've solved it!
Working code is:
a0 = [0 0 0];
b0 = a0/3 .* [1 1 1];
a1 = a0/2 .* [0 1 1];
a2 = a0/2 .* [1 0 1];
a3 = a0/2 .* [1 1 0];
In = n .*a1 + m .*a2 + l .*a3;
P = n .*a1 + n .*a2 + l .*a3 + b0;
for n = 1:5
for m = 1:5
for l = 1:5
In(n+1,m+1,l+1,1:3) = a0';
P(n+1,m+1,l+1,1:3) = a0' + b0';
end
end
end