function [x,w]=gaussQuad(n)
% Computes interpolation points x_j and weights w_j of the Gauss quadrature

b = 1:n; b = b./sqrt(4*b.*b-1);
J=diag(b,-1)+diag(b,1); [ev,ew]=eig(J);
x=diag(ew); w=(2*(ev(1,:).*ev(1,:)));
