# File "Rnewt".

g = function(theta) {+ 4*(1-theta)^3 + 4*(2-theta)^3 }

gp = function(theta) {- 12*(1-theta)^2 - 12*(2-theta)^2 }

xlist = c(4)

for (i in 1:100) {

    xprev = xlist[length(xlist)]

    xlist = c(xlist, xprev - g(xprev) / gp(xprev))

}

plot(xlist)