A regression with interaction terms can be powerful but it can also be difficult to interpret.

For one thing, at what level do you expect interactions? You can tell R to only look at first-order interactions using a formula such as y~(x1+x2+x3)^2, or you can tell it to look at all possible interactions with y~ x1*x2*x3, or even limit interactions using y~x1*(x2+x3). Using too many risks overfitting.

Then you have the model selection problem - how do you identify and prune out interactions you don't need? I'd use the 'drop1' command in R, which will give you the AIC resulting from dropping terms, and it emphasizes dropping interactions before dropping main effects.

As all of these terms you listed are continuous, you will want to center them by subtracting means or medians. That will make coefficients easier to interpret. Can we assume linearity? Perhaps your most appropriate model isn't y~x1*x2 but y~ln(x1)*x2^2 for example.

I realize this is a lot. As a first step, graph. If you converted two of the three X's you list into factors, say with 2 to 4 levels, then you can plot Y vs (continuous) X1, grouped by (factor) X2 within each plot, and using subplots based on (factor) X3. This is easy to set up in ggplot2. If there are meaningful interactions or non-linearities, they should show up. You may also want to explore contour graphs.

Other potentially useful tools:

Regression trees such as 'rpart'. The branching patterns may indicate interactions present in the data.

The 'party' package combines regression trees with regression models. The 'mob' command in that package, as I recall, generates a regression tree terminating in linear regression models.

Good luck!