I'm running the following code on quite large data frames. I've rewritten it for the iris dataset to make it reproducable.
I'm quite unexperienced with the apply functions and I find them a pain in the bum to apply them.

Is there any ways to drastically improve the performance of this process?

lmfit <- lm(iris$Petal.Width ~ iris$Sepal.Length + iris$Sepal.Width)
out_index <- 1
TableWithResiduals <- data.frame(matrix(ncol = ncol(iris) +1, nrow = nrow(iris)))
for (row in 1:length(resid(lmfit))){
  TableWithResiduals[out_index,] <- cbind(iris[row,],resid(lmfit)[row])
  out_index <- out_index +1  
}
colnames(TableWithResiduals) <- colnames(iris)
colnames(TableWithResiduals)[length(TableWithResiduals)] <- "Residual_value"