x <- cbind(seq(1, 50, 1), seq(51, 100, 1))

y <- x[,1]*x[,2]

colnames(x) <- c('x1', 'x2')
names(y) <- 'y'

dt <- data.frame(x, y)

model <- neuralnet(y~., dt, hidden=10, threshold=0.01)

n <- names(dt)
f <- as.formula(paste("y ~", paste(n[!n %in% "y"], collapse = " + ")))
f

## gives
> f
y ~ x1 + x2

## fit model using `f`
model <- neuralnet(f, data = dt, hidden=10, threshold=0.01)

> model
Call: neuralnet(formula = f, data = dt, hidden = 10, threshold = 0.01)

1 repetition was calculated.

        Error Reached Threshold Steps
1 53975276.25     0.00857558698  1967

f <- reformulate(setdiff(colnames(dt), "y"), response="y")