function [ptwt] = legptwt(N); format long e xn=[1:1:N]; xnsq=xn.*xn; % beta_n recurrence coefficients b=xnsq./(4*xnsq-1); rtb=sqrt(b); rtb(N)=[]; % Jacobi matrix t=diag(rtb,-1)+diag(rtb,1); [f,lambda]=eig(t); % quadrature points as the eigenvalues pt=diag(lambda); % quadrature weights from the 1st component of the ith eigenvector wt=2*f(1,:).^2; display('Legendre quadrature points and weights') for n=1:N fprintf('%2i %26.16f %26.16e\n', n, pt(n),wt(n)) end %ptwt=[pt,wt'];