function q = composite_newton_cotes ( m, n, a, b, f )

  breaks = linspace ( a, b, m + 1 );

  q = 0.0;
  for i = 1 : m
    ai = breaks(i);
    bi = breaks(i+1);
    q = q + newton_cotes ( n, ai, bi, f );
  end

  return
end
