This function is typically used inside a model formula
  in order to specify the most general possible polynomial
  of order n involving one numerical variable x
  or two numerical variables x,y.
It is equivalent to poly(, raw=TRUE).
If only one numerical vector argument x is given, 
  the function computes the vectors x^k for
  k = 1, 2, …, n. These vectors are combined into a matrix
  with n columns.
  
If two numerical vector arguments x,y are given,
  the function computes the vectors x^k * y^m for
  k >= 0 and m >= 0 satisfying
  0 < k + m <= n. These vectors are combined into a matrix
  with one column for each homogeneous term.