R function: solve

m   <- c(5, 3, 2, 1,
         7,10,-1, 4,
         4, 6, 6, 7,
        -5,-2, 0, 6)

dim(m) <- c(4, 4)

n <- solve(m)

round(n %*% m, 15)
#      [,1] [,2] [,3] [,4]
# [1,]    1    0    0    0
# [2,]    0    1    0    0
# [3,]    0    0    1    0
# [4,]    0    0    0    1