Generating Polynomial Features in R
You can use two methods: one that takes a formula argument or if you are not using formula you can use model.matrix.

You can do the following:

# Using the same data from the question:
data &lt;- data.table(A = c(0, 2, 4), 
                   B = c(1, 3, 5))

&gt; data

   A B
1: 0 1
2: 2 3
3: 4 5


If only the raw data and the interactions are required

formula = y ~ .^2 

model.matrix(formula, data=data)
  (Intercept) A B A:B
1           1 0 1   0
2           1 2 3   6
3           1 4 5  20


If raw data, squared columns, interactions and intercept are needed, just like the question

You can use the following formula for all variables:

formula = y ~ .^2 + poly(var, 2, raw=TRUE)[, 2] ... etc


In addition you can paste the variables names automatically. Formula based on this post

formula &lt;- as.formula(paste(' ~ .^2 + ',paste('poly(',colnames(data),',2, raw=TRUE)[, 2]',collapse = ' + ')))

&gt; formula 

~.^2 + poly(A, 2, raw = TRUE)[, 2] + poly(B, 2, raw = TRUE)[, 2]

&gt; model.matrix(formula, data=data)

  (Intercept) A B poly(A, 2, raw = TRUE)[, 2] poly(B, 2, raw = TRUE)[, 2] A:B
1           1 0 1                           0                           1   0
2           1 2 3                           4                           9   6
3           1 4 5                          16                          25  20


The resulting dataframe is identical to the array in the question.

If squared columns, interactions and intercept are needed.

&gt; formula &lt;- as.formula(paste(' ~ A:B + ',paste('poly(',colnames(data),',2, raw=TRUE)[, 2]',collapse = ' + ')))

&gt; formula
~A:B + poly(A, 2, raw = TRUE)[, 2] + poly(B, 2, raw = TRUE)[, 2]

&gt;  model.matrix(formula, data=data)
  (Intercept) poly(A, 2, raw = TRUE)[, 2] poly(B, 2, raw = TRUE)[, 2] A:B
1           1                           0                           1   0
2           1                           4                           9   6
3           1                          16                          25  20