The error message you see "Function '[' is not in the derivatives table" is because D can only recognize a certain set of functions for symbolic operations. You can find them in ?D:
The internal code knows about the arithmetic operators ‘+’, ‘-’,
‘*’, ‘/’ and ‘^’, and the single-variable functions ‘exp’, ‘log’,
‘sin’, ‘cos’, ‘tan’, ‘sinh’, ‘cosh’, ‘sqrt’, ‘pnorm’, ‘dnorm’,
‘asin’, ‘acos’, ‘atan’, ‘gamma’, ‘lgamma’, ‘digamma’ and
‘trigamma’, as well as ‘psigamma’ for one or two arguments (but
derivative only with respect to the first). (Note that only the
standard normal distribution is considered.)
While the "[" is actually a function in R (read ?Extract or ?"[").
To demonstrate the similar behaviour, consider:
s <- function (x) x
D(expression(s(x) + x ^ 2), name = "x")
# Error in D(expression(s(x) + x^2), name = "x") :
# Function 's' is not in the derivatives table
Here, even though s has been defined as a simple function, D can do nothing with it.
Your problem has been solved by my recent answers for Function for derivatives of polynomials of arbitrary order (symbolic method preferred). Three methods are provided in three of my answers, none of which are based on numerical derivatives. I personally prefer to the one using outer (the only answer with LaTeX math formula), as for polynomials everything is exact.
To use that solution, use the function g there, and specify argument x by values where you want to evaluate the derivative (say 0:10), and pc by your polynomial regression coefficients s. By default, nderiv = 0L so the polynomial itself is returned as if predict.lm(m1, newdata = list(a = 0:10)) were called. But once nderiv is specified, you get exact derivatives of your regression curve.
a <- 0:10
b <- c(2, 4, 5, 8, 9, 12, 15, 16, 18, 19, 20)
plot(a, b)
m1 <- lm(b ~ a + I(a ^ 2) + I(a ^ 3))
s <- coef(m1)
#(Intercept) a I(a^2) I(a^3)
# 2.16083916 1.17055167 0.26223776 -0.02020202
## first derivative at your data points
g(0:10, s, nderiv = 1)
# [1] 1.1705517 1.6344211 1.9770785 2.1985237 2.2987568 2.2777778 2.1355866
# [8] 1.8721834 1.4875680 0.9817405 0.3547009
Other remark: It is suggested that you use poly(a, degree = 3, raw = TRUE) rather than I(). They do the same here, but poly is more concise, and make it easier if you want interaction, like in How to write interactions in regressions in R?
