I have written a code but I am stuck at a point after finding the tensor product. I need to find basis of chebyshev polynomial after tensor product . using: shapefcn = polyBasis(​'chebyshev​',degree) would provide the correct basis functions ?

m = 3; % chebyshev Polynomial of order m+1 n = 6; % Number of Random Variables

total_exp0 = (m+1)^n;

p = factorial(n+m)/(factorial(n)*factorial(m))-1;

pcoeff = ChebT(m+1); %genertaes Chebyshev polynomial coefficients zeta_c = roots(pcoeff);

%Tensor product of the roots

quad_nodes0 = zeros((m+1)^n,n);

for c=1:(m+1)^n

    temp=dec2base(c-1,m+1,n);

    for d=1:n

        quad_nodes0(c,d)=zeta_c(str2double(temp(d))+1);

    end

end

if ((p+1)*10 < total_exp0)

   total_exp=(p+1)*10;

   tmprows=randperm(total_exp0,total_exp);

   % tmprows=I(length(I)-total_exp+1:length(I));

   quad_nodes=quad_nodes0(tmprows,:);

else

   total_exp=total_exp0; 

   quad_nodes=quad_nodes0;

end

ExpMat = zeros(total_exp,p+1);