Confidence interval for GLM or the maximum of a function? Imagine I have a set of (xi,yi) measures.
I can show it on a scatter plot

I want to choose the value of x that maximizes  y,
or I could fit a function and find the values of the parameters that maximize that function.
To put it simple I decide to fit a general linear regression model  (quadratic on the variable) 
Y=a+bx+cx^2.
I can do it with R, using the glm()  funtcion.
And I would get the values of a, b, c their standard errors  SE(a), SE(b), SE(c), R2 and more things.
And the max would be located at x=-b/2c  
My question is,  How can I calculate the confidence interval for that x, the predictor?
Is it just the typical formula   x+-Z·S/sqrt(n)?
and using some result provided by R for that S?
or do I need to calculate that standard using a more complex methods, such as the square root of the sum of derivatives squared...?
Or it's something completly different?
How would you do it with a simulation with x ?
maybe calculate "y" for different values of x with perturbation. 
then choose the 95% max values of y (??) and calculate the confidence interval for the associated x's ?
How would you do it?
 A: In the comments there are good proposals, like the traditional delta method. A newer method which might be more exact, is constructing a profile likelihood confidence interval, which today is quite easy.
I will show below with an simulated example, in R. First write the linear regression (polynomial) model as
$$ Y_i = a + b x_i + c x_i^2 + \epsilon_i $$
Now, assuming $c<0$ so as to have a maximum ... (If $c>0$ the same methods will give a confidence interval for the minimum), the maximum is at
$$ M=-\frac{b}{2c} $$
$M$ is our interest (or focus) parameter, so reparametrize the regression model as
$$ Y_i = a + b x_i\left( 1-\frac{x_i}{2M}\right) +\epsilon_i $$
This is now no longer linear in the parameters, but we can fit it as a nonlinear regression model:
a <- 1
b <- 1
c <- -1/2
x <- seq(from=-5,  to=5,  by=1/3)

set.seed(7*11*13) # My public seed
Y <- a + b*x  +  c*x^2  +  rnorm(length(x), 0, 2)

mydata <- data.frame(Y, x)

mod0 <- nls(Y ~ a  + b*x*(1-x/(2*M)), data=mydata,
            start=list(a=0, b=0.5, M=2))  

 summary(mod0)

Formula: Y ~ a + b * x * (1 - x/(2 * M))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a   1.2778     0.6284   2.033   0.0516 .  
b   0.9671     0.1404   6.888 1.74e-07 ***
M   0.8975     0.1571   5.711 3.99e-06 ***
  ...  

Now we call the confint function, which calls the profile function , doing likelihood profiling, and constructs the confidence interval from there:
 confint(mod0, 3)
Waiting for profiling to be done...
     2.5%     97.5% 
0.6044179 1.2658407 

We can also plot the likelihood profile for M:

Now, as an exercise readers can compare this with the delta method ...
