function q = simpson13 ( a, b, f )

  x1 = a;
  w1 = ( b - a ) * 1.0 / 6.0;

  x2 = ( a + b ) / 2.0;
  w2 = ( b - a ) * 4.0 / 6.0;

  x3 = b;
  w3 = ( b - a ) * 1.0 / 6.0;

  q = w1 * f ( x1 ) + w2 * f ( x2 ) + w3 * f ( x3 );

  return
end
