Creating a function called factorize which takes a square matrix and decomposes into L and U matrices.

factorize <- function(A) 
{
  if (nrow(A) == ncol(A)) 
  {
      size <- nrow(A)
      # Constructing a diagonal matrix. This will become the matrix L after factorization.
      L <- diag(size) 
      for (i in 1:(size - 1)) 
      {
          for (j in (i + 1):size) 
          {
              L[j, i] <- A[j, i] / A[i, i]  # building the L matrix
              A[j, ]  <- A[j, ] - L[j, i] * A[i, ] # building U matrix
          }
      }
      
      print(L %*% A)
      
      print('-----L--------')
      print(L)
      print('--------U--------')
      print(A)
      
  } 
}

testing with a 3 X 3 Matrix.

A <- matrix(c(3,5,6,4,5,1,5,1,4), nrow = 3)
print(A)
##      [,1] [,2] [,3]
## [1,]    3    4    5
## [2,]    5    5    1
## [3,]    6    1    4
factorize(A)
##      [,1] [,2] [,3]
## [1,]    3    4    5
## [2,]    5    5    1
## [3,]    6    1    4
## [1] "-----L--------"
##          [,1] [,2] [,3]
## [1,] 1.000000  0.0    0
## [2,] 1.666667  1.0    0
## [3,] 2.000000  4.2    1
## [1] "--------U--------"
##      [,1]      [,2]      [,3]
## [1,]    3  4.000000  5.000000
## [2,]    0 -1.666667 -7.333333
## [3,]    0  0.000000 24.800000