function poly_gauss ( )

%% poly_gauss uses Gauss quadrature to estimate integral poly(x) over [-1,3].
%
  a = -1.0;
  b = 3.0;
  q_exact = poly_int ( a, b );

  for n = 1 : 5
    [ xi, rho ] = gauss ( n );
%
%  Transform the rule from [-1,+1] to [a,b].
%
    x = a + ( xi + 1.0 ) * ( b - a ) / 2.0;
    w = ( b - a ) * rho / 2.0;

    q = w' * poly ( x );
    fprintf ( 1, '  %d  %g  %g\n', n, q, q - q_exact );
  end

  return
end

function y = poly ( x )
  y = 10*x.^4 - 12*x.^3 + 2*x;
  return
end

function value = poly_int ( a, b )
  ya = 2*a.^5 - 3*a.^4 + a.^2;
  yb = 2*b.^5 - 3*b.^4 + b.^2;
  value = yb - ya;
  return
end

