function q = simpson38 ( a, b, f )

  x1 = ( 3.0 * a + 0.0 * b ) / 3.0;
  w1 = ( b - a ) * 3.0 / 24.0;

  x2 = ( 2.0 * a + 1.0 * b ) / 3.0;
  w2 = ( b - a ) * 9.0 / 24.0;

  x3 = ( 1.0 * a + 2.0 * b ) / 3.0;
  w3 = ( b - a ) * 9.0 / 24.0;

  x4 = ( 0.0 * a + 3.0 * b ) / 3.0;
  w4 = ( b - a ) * 3.0 / 24.0;

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

  return
end
