> help d_mid

  uses the midpoint rule to estimate f'(x) 
 
  x is the point the derivative is being approximated at, h is the 
  interval distance the approximation is over.  Requires suboroutine
  f.m
 
  function y = d_mid(x,h)

> x=3;
> y_true = -sin(x);

approximate the derivative of cos(x) at x=3, using intervals that
get smaller by factors of 2. 

> y1 = d_mid(x,1/2);
> y2 = d_mid(x,1/4);
> y3 = d_mid(x,1/8);
> y4 = d_mid(x,1/16);
> y5 = d_mid(x,1/32);
> format short e
> format compact
> y_true - y1
ans =
 -5.8069e-003
> y_true - y2
ans =
 -1.4654e-003
> y_true - y3
ans =
 -3.6721e-004
> y_true - y4
ans =
 -9.1857e-005
> y_true - y5
ans =
 -2.2968e-005

Since we know the derivative of cos(x), we know the error and we can
test the convergence rate directly through the ratios of the errors:
> (y_true - y1)/(y_true-y2)
ans =
  3.9627e+000
> (y_true - y2)/(y_true-y3)
ans =
  3.9906e+000
> (y_true - y3)/(y_true-y4)
ans =
  3.9977e+000
> (y_true - y4)/(y_true-y5)
ans =
  3.9994e+000
The above is consistent with an error that scales like h^2

Now we repeat the above with intervals that get smaller by factors
of 3.
> y1 = d_mid(x,1/2);
> y2 = d_mid(x,1/(2*3));
> y3 = d_mid(x,1/(2*9));
> y4 = d_mid(x,1/(2*27));
> y5 = d_mid(x,1/(2*81));

This time, the ratios of the errors go to 9, again consistent with
an error that scales like h^2
> (y_true - y1)/(y_true-y2)
ans =
  8.9005e+000
> (y_true - y2)/(y_true-y3)
ans =
  8.9889e+000
> (y_true - y3)/(y_true-y4)
ans =
  8.9988e+000
> (y_true - y4)/(y_true-y5)
ans =
  8.9999e+000


If you don't know the derivative then you don't know the error.  In
general, we don't know the derivative (why else would we be asking
matlab to help us find it?)  But we can still use the approximatations
to estimate the rate that error decreases to zero.
> (y1-y2)/(y2-y3)
ans =
  8.8895e+000
> (y2-y3)/(y3-y4)
ans =
  8.9877e+000
> (y3-y4)/(y4-y5)
ans =
  8.9986e+000
> diary off