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.