function [x, c] = gaussquad(n)
% Compute Gaussian quadrature points and coefficients.
% For education only: the monomial basis is unstable for high n.

P = zeros(n + 1, n + 1);
P([1, 2], 1) = 1;
for k = 1:n - 1
    P(k + 2, 1:k + 2) = ((2 * k + 1) * [P(k + 1, 1:k + 1), 0] - ...
        k * [0, 0, P(k, 1:k)]) / (k + 1);
end
x = sort(roots(P(n + 1, 1:n + 1)));

A = zeros(n, n);
for i = 1:n
    A(i, :) = polyval(P(i, 1:i), x)';
end
c = A \ [2; zeros(n - 1, 1)];
