function [xg,wg] = findgauss(n)
%    [xg,wg] = findgauss(n)
%    find Gauss nodes xg(1),...xg(n) and weights wg(1),...,wg(n) for [-1,1]
%
%  NOTE: findgauss illustrates the method of finding Gauss nodes and weights
%        explained in class. For computations you should use [xg,wg] = gauleg(n) 
%        instead. This is more accurate and efficient.
%
%    The integral of f(x) over x from -1 to 1
%    can be approximated by   
%      Q = wg(1)*f(xg(1)) + ... + wg(n)*f(xg(n))
%
%    The integral of f(x) over x from a to b
%    can be approximated by   
%
%      tg = (a+b)/2 + xg*(b-a)/2;
%      Q = (b-a)/2*( wg(1)*f(tg(1)) + ... + wg(n)*f(tg(n))
%
% Example: Approximate integral of cos(x) for x from 1 to 4 
%          using Gaussian integration with 3 nodes
%
%      [xg,wg] = findgauss(3)
%      a = 1;  b = 4;
%      tg = (a+b)/2 + xg*(b-a)/2;
%      Q = (b-a)/2 * wg'*cos(tg)
 
% determine ex(k+1) = integral of x^k from -1 to 1 for k = 0,...,2*n-1
for k=0:2*n-1
  ex(k+1) = ( (1)^(k+1) - (-1)^(k+1) ) / (k+1);
end

% STEP 1: Find polynomial omega(x) = x^n + c(n)*x^{n-1} + ... + c(1)
%         such that integral of omega(x)*x^k from -1 to 1 is 0 for k=0,...,n-1
%         gives linear system A*c = b
for j=1:n
  A(j,:) = ex(j:n+j-1);
end
b = - ex(n+1:2*n)';
c = A\b;

% STEP 2: Gauss nodes xg(1),...,xg(n) are roots of polynomial omega(x) 
%         which has coefficients [ 1; c(n); c(n-1); ... ; c(1) ] = [ 1; c(end:-1:1) ]              

xg = roots([ 1; c(end:-1:1) ]);

% STEP 3: Determine weights wg(1),...,wg(n) 
%         such that Gauss rule is exact for f(x)=1,...,x^{n-1}
%         gives linear system M*wg = y
for j=1:n
  M(j,:) = xg'.^(j-1);
end
y = ex(1:n)';
wg = M\y;      
