When I change the polynomial to:
cf = fliplr(coeffs(x*3 + 2*z, [x y z],'All'))
I get a 3-D array for the coefficients, as expected. What you see might be the effect of matlabs reduction of the trailing dimensions. Compare for example with:
size(randn(12,23,1))
I don't know how to solve this for you, unfortunately. Maybe you can add a term with high enough degree in your last variable to avoid that type of reduction and then crop the output for your future manipulations - this seems like a very messy idea.
HTH
