how to add second order terms into the model in R? My data mydata consists of columns of x1,x2,..,x100,y in R. But I am thinking a linear model with second order terms such as y ~ x1^2 + x2^2 + x1*x2 + ... how do I achieve that within formula or in any other way in R?
When I tried above, my model pls ignored all second order terms. Do I have to manually create those columns?
 A: The formula documentation for R shows how to do this.  In short, you use poly().  For example, make some quadratic data:
x <- rnorm(100)
y <- x + x**2 * 0.5 + rnorm(100)  

Now fit this using a second order polynomial (i.e x and x**2) like this
mod <- lm(y ~ poly(x, 2))

Note that this will fit an orthogonal polynomial, so it won't recover 1 and 0.5 as the coefficients in the generating distribution.  If for some reason you want that, use poly(x, 2, raw=TRUE).  In general you don't for stability reasons, so stick with the cooked version.
There is also polym as in: lm(y ~ polym(x, z, degree=2) for a model with a full set of crossed variables, which is a bit more trouble to interpret, but that's presumably not important with hundreds of variables.
A: Type : 
lm(y ~ x1 + x2 + I(x1*x2) + I(x1^2) + ...)

A: Here's how to do it in principle, illustrated on a smaller dataset with only 10 predictors:
# Make fake data
mydata = as.data.frame(matrix(rnorm(1100), 100))
names(mydata) = c(paste0("x", 1:10), "y")

# Form a matrix containing all predictor columns but not y
x = as.matrix(mydata[, 1:10])
# Use poly() to form all 2-way interactions and 2nd order terms
x2 = poly(x, degree = 2, raw = TRUE)
# Resave as a data frame including y
mydata2 = as.data.frame(cbind(x2, y = mydata$y))
# Fit the complete linear model
lm2 = lm(y ~ ., data = mydata2)

However, you have 100 predictors. In my experience, with more than 10-15 predictors, R usually cannot allocate enough memory for the matrix containing every 2-way interaction. You will get unhelpful errors or R will simply crash.
If so, consider whether you really need all 2-way interactions. Maybe just a subset would make sense. For instance, you could use poly() as above to form all 2-way interactions within one subset of x's, then again to form interactions between another subset of x's, but not have any interactions across those subsets.
