source('../Chapter3/utils.R')

## The eigenvalues are the roots of the following function:
##
eigenvalue_polynomial = function(x){
    x^4 - 6*x^2 - 8*x - 3
}
polynomial_root_bounds(c(-3, -8, -6, 0, 1))

## Lets get a feel for where the zeros of this polynomial i.e.:
##
xs = seq(-3, +4, length.out=1000)
ys = sapply(xs, eigenvalue_polynomial)
plot(xs, ys, 'l', ylab='p(x)')
grid()

## Compute the roots with Muller's method:
##
x0s = c(-3, 0, 3.5)
res = mullers_iteration(x0s, eigenvalue_polynomial, 4, xtol=1.e-6, ftol=1.e-3, verbose=FALSE)
print(res$x)

## Verify our results using polyroot:
##
print(polyroot(c(-3, -8, -6, 0, 1)))