How to construct Legendre polynomials for $x_1,...,x_k$?

Question

I'm trying to run a nonparametric regression to estimate the unknown conditional mean $E(Y|X_1=x^*_1,X_2=x^*_2)$ using data set $\{Y_i,X_{1i},X_{2i}\}_{i=1}^n$. This could be done by nonparametric series regression: first construct $K$-th order polynomial regressors $r_1(x_1,x_2),...,r_K(x_1,x_2)$ and then regress $Y_i$ on $r_1(X_{1i},X_{2i}),...,r_K(X_{1i},X_{2i})$ and use the predicted value at $(x_1^*,x_2^*)$ as our estimate. I plan to use the $K-$th order Legendre polynomials in $x_1,x_2$,and I'm wondering how shoud I construct them? From linear algebra notes I found that if I only have one variable $x$ then the up-to-3rd order Legendre polynomials are: 1, $x$, $\frac{3x^2-1}{2}$, $\frac{5x^3-3x}{2}$. I'm wondering what are the corresponding Legendre polynomials when I have $x_1,x_2$? A formula, recursive formula, reference or software package/function that generates this are all ok. Thanks