Newton Raphson in R

> library(polynom)
> coefs <- c(-1, -2, 3, -3, -2, 1)
> p <- polynomial(coef = coefs)
> p1 <- deriv(p)
> p.f <- as.function(p)
> p1.f <- as.function(p1)
> x <- 0
> y <- p.f(x)
> xs <- vector()
> for (i in 1:6) {
+     x <- x - y/p1.f(x)
+     y <- p.f(x)
+     print(paste(x, y, sep = "   "))
+     xs <- c(xs, x)
+ }
[1] "-0.5   0.96875"
[1] "-0.336842105263158   0.0986450190238757"
[1] "-0.315728448396289   0.00191853843116618"
[1] "-0.315301162703277   7.90693466257153e-07"
[1] "-0.315300986459363   1.34559030584569e-13"
[1] "-0.315300986459333   0"

> par(mfrow = c(1, 1))
> plot(p)
> abline(h = 0)
> abline(v = 0)

> print(xs, digits = 16)
[1] -0.5000000000000000 -0.3368421052631579 -0.3157284483962894
[4] -0.3153011627032768 -0.3153009864593627 -0.3153009864593327
> solve(p)[2]
[1] -0.3153009864593327+0i