source('utils.R')

##
## Assume we have estimates of the eigenvalues from running the R code section_4_8_ex_13.R
##
evectors = c()
for( ev in newton_roots ){
    print(sprintf('Computing the eigenvector of B corresponding to the eigenvalue= %10.6f', ev))
    B_minus_pI = B - (ev+1.e-2)*diag(n)
    B_minus_pI_inv = solve(B_minus_pI, diag(n)) ## the "power method" is applied to this matrix
    res = power_method(B_minus_pI_inv) 
    evectors = c(evectors, res$max_eigenvector)
}
EV = matrix(data=evectors, ncol=length(newton_roots), nrow=length(newton_roots), byrow=FALSE)
print('Eigenvalues are= ')
print(EV)