I went into the Simpsons program and introduced a programming error
by making the for loop stop too early.  Below, I show that this makes
the program no better than the left-hand rule.

> [a,h,e1] = simp(1,2,8+1);
> [a,h,e2] = simp(1,2,16+1);
> [a,h,e3] = simp(1,2,32+1);
> format short e, format compact
> e1
e1 =
 -7.4689e-002
> e2
> e2 =
 -4.4788e-002
> e3
> e3 =
 -2.4217e-002
> e1/e2
ans =
  1.6676e+000
> e2/e3
ans =
  1.8495e+000
> [a,h,e4] = simp(1,2,64+1);
> [a,h,e5] = simp(1,2,128+1);
> e3/e4
ans =
  1.9283e+000
> e4/e5
ans =
  1.9650e+000

Now I fixed my Simpsons program so it wasn't buggy and redid the 
computation.  This time, I computed the rate of convergence from
the approximate areas instead of from the errors.  I get 16, as
desired.

> [a1,h,e1] = simp(1,2,8+1);
> [a2,h,e2] = simp(1,2,16+1);
> [a3,h,e3] = simp(1,2,32+1);
> [a4,h,e4] = simp(1,2,64+1);
> [a5,h,e5] = simp(1,2,128+1);

> (a1-a2)/(a2-a3)
ans =
    1.602347475658662e+001
> (a2-a3)/(a3-a4)
ans =
    1.600586622522044e+001
> (a3-a4)/(a4-a5)
ans =
    1.600138798725372e+001


Now we move on to differentiation.  We want to approximate the
derivative of cos(x) at x=2.

> x = 2;
> help d_right

  uses the right-hand 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_right(x,h)

> y1 = d_right(x,1);
> y2 = d_right(x,1/2);
> y3 = d_right(x,1/4);
> y4 = d_right(x,1/8);
> y1
y1 =
   -5.738456600533031e-001
> y2
y2 =
   -7.699935579995826e-001
> y3
y3 =
   -8.481071447023869e-001
> y4
y4 =
   -8.809559852573012e-001
> format short e

We compare the approximate derivative to the true derivative:
> y_true = -sin(x)
> y_true =
 -9.0930e-001
> y_true-y1
ans =
 -3.3545e-001
> y_true-y2
ans =
 -1.3930e-001
> y_true-y3
ans =
 -6.1190e-002
> y_true-y4
ans =
 -2.8341e-002

We compute the ratio of the errors, and find that they are close
to 2, as expected.
> (y_true-y1)/(y_true-y2)
ans =
  2.4081e+000

Now we repeat the exercise using the mid-point rule for the derivative
and find that the ratio looks closer to 4...
> y1 = d_mid(x,1);
> y2 = d_mid(x,1/2);
> y_true-y1
ans =
 -1.4415e-001
> y_true-y2
ans =
 -3.7417e-002
> (y_true-y1)/(y_true-y2)
ans =
  3.8526e+000
> diary off