When we want to create a model with an arbitrary number of independent variables, we can use the below:

create_lm <- function(data, dep, covs) {
# Create the first part of the formula with the dependent variable
form_base <- paste(dep, "~")
# Create a string that concatenates your covs vector with a "+" between each variable
form_vars <- paste(covs, collapse = " + ")
# Paste the two parts together
formula <- paste(form_base, form_vars)
# Call the lm function on your formula
lm(formula, data = data)
}

For example, using the built-in mtcars dataset:

create_lm(mtcars, "mpg", c("wt", "cyl"))
Call:
lm(formula = formula, data = data)
Coefficients:
(Intercept) wt cyl
39.686 -3.191 -1.508

The downside is that the printed output from the model doesn't reflect the particular call you made to lm, not sure if there is any way around this.