How to build an 2-D polynomial from 1-D orthogonal polynomials

Question

I have an set of orthogonal polynomials such as

I want to build an 2D polynomial following the equation $$P_k(x,y)=P_k(x)P_k(y)$$ where $k=1..4, (x,y) \in [-1, 1]^2$

Based on given $P_n(x)$ as above. How to construct an 2-D $P_k(x,y)$? I look at one solution but I did not clear understand. If you understand it, please explain it help me

%% ten 2-D basic function
bais(:,:,1) = 1;
bais(:,:,2) = x;
bais(:,:,3) = (3.*x.*x - 1)./2;
bais(:,:,4) = (5.*x.*x.*x - 3.*x)./2;
bais(:,:,5) = y;
bais(:,:,6) = x.*y;
bais(:,:,7) = y.*(3.*x.*x -1)./2;
bais(:,:,8) = (3.*y.*y -1)./2;
bais(:,:,9) = (3.*y.*y -1).*x./2;
bais(:,:,10) = (5.*y.*y.*y -3.*y)./2;
B = bais; 
for kk=1:10
    A=bais(:,:,kk).^2;
    r = sqrt(sum(A(:)));
    B(:,:,kk)=bais(:,:,kk)/r;
end