how to compute maximum of fitted polynomial?

Thank you all!  This approach, using the 'polynom' library, did the trick.

library(polynom) # -6 + 11*x - 6*x^2 + x^3
p0 <- polynomial(c(-6, 11, -6, 1)) # has zeros at 1, 2, and 3
p1 <- deriv(p0); p2 <- deriv(p1) # first and second derivative
xm <- solve(p1) # maxima and minima of p0
xmax = xm[predict(p2, xm) < 0] # select the maxima
xmax # [1] 1.42265

With these tweaks:

In fitting the model to the data I had to use raw=TRUE:
model <- lm( y ~ poly(x, 3, raw=TRUE) )
 
Then I could generate p0 directly from my lm:
p0 <- polynomial( coef(model) )

And I answered my second question by using the obvious:
predict(p2,xmax)

I don't know what I would have done if the optima weren't between x=0 and x=1, which was my constraint.  In that case the "maximum" would have been one of the endpoints rather than a zero of p1.  I suppose I could just have checked for it with some if/then code.  Fortunately it didn't turn out to be an issue with my data.

And no, this wasn't a homework problem.  I didn't do the math by hand because I needed to automate this process for several subsets of my data and several fitted models.