  n = 7
  [ xi, rho ] = gauss ( n )
  xi_sum = sum ( xi )
  rho_sum = sum ( rho )
  xi_sym = xi(1:n) + xi(n:-1:1)
  rho_sym = rho(1:n) - rho(n:-1:1)

  f = @(x)cos(x);
  f_anti = @(x)sin(x);
  gn = rho' * f(xi)
  in = f_anti(+1.0) - f_anti(-1.0)

  f = @(x) 13*x.^12-6*x.^5+5*x.^4-7;
  f_anti = @(x) x.^13-x.^6+x.^5-7*x;
  gn = rho' * f(xi)
  in = f_anti(+1.0) - f_anti(-1.0)

