lm() and interactions in model formula for x passed as matrix

Suppose I have x variables x1, x2, x3 (however in general I don't know
how many x variables there are). I can do
X<-cbind(x1,x2,x3)
lm(y ~ X)
This fits the no-interaction model with b0, b1, b2, b3.

How can I get lm() to fit the model that includes interactions when I
pass X to lm()? For my example,
lm(y~x1*x2*x3)
I am looking for something along the lines of
lm(y~X ...)
where ... is some extra stuff I need to fill in.

The formula syntax in R allows you to specify interactions with the  
"^" operator but some testing makes me think you cannot use either y  
~ .^3 or y ~ X^3 with matrix data arguments, here assuming you only  
want interaction up to third order.

Assuming you know how to use do.call("cbind", varlist)
perhaps:

  form = as.formula( paste("y ~ (",
                            paste(colnames(X), collapse="+"),
                              ")^3", sep="")  )
lm(form)


--- output------
Call:
lm(formula = form)

Coefficients:
(Intercept)           x1           x2           x3        x1:x2         
x1:x3
   -0.383296    -0.333429     0.003976     0.332982    -0.001130      
0.100698
       x2:x3     x1:x2:x3
    0.366745     0.122111

David Winsemius, MD
West Hartford, CT