Example

> f:=x->ln(cos(x));

> df:=D(f):ddf:=D(df):ddf(0.4);

> for i from 1 to 3 do h:=10.^(-i): ddf1:=(f(0.4-h)-2*f(0.4)+f(0.4+h))/h^2: acterr:=abs(ddf1-ddf(0.4)): print(h,ddf1,acterr, acterr/h^2); od:

This shows that the formula does work, but care must be taken for round-off errors. Again, an optimal value of will exist .