"Toeplitz" <- function(x) { matrix(x[1 + abs(outer(seq(along = x), seq(along = x), FUN = "-"))], byrow = T, ncol = length(x)) } "test" <- function(phi,n){ solve(Toeplitz(phi^seq(0,n-1))) } print(system.time(test(0.5, 1000))[1])