Not sure what you mean...

in R you can use ~(A+B+C+..)^k to model all k-way interactions of the predictors A, B, C...

Example:

Y ~ (A+B+C+D)^3

is the model including all 2-way and 3-way interactions of the 4 predictors:

Y ~ A+B+C+D+A:B+A:C+A:D+B:C+B:D+C:D+A:B:C+A:B:D+A:C:D+B:C:D

If you have calculated a model with all these interactions you can save the output from the summary-method in a variable and filter the interaction-coefficients by searching for the colon (":") in the row names of the coefficients table.

I understand that in R is

response ~ A*B*C*...*N

It produces a model with each single variable, each pair, and all remaining interactions. Of course, it is needed a lot of sample size for many variables.

If you want to fit a model with all interactions automatically included, the question has been answered abundantly above. In case what you want is to select relevant interaction terms (" check all combination of interactions"), then you might want to use something like function stepAIC in R package MASS.  Assuming your variables are y, x1, x2 and x3 defined as columns in the data frame mydata, a command such as:

model