function [alpha,beta,b,P]=OrthogonalPolynomialCoefficients(x,y)
% ORTHOGONALPOLYNOMIALCOEFFICIENTS 3 term recurrence coefficients
%   [alpha,beta,b,P]=OrthogonalPolynomialCoefficients(x,y); computes
%   the coefficients of the 3-term recurrence of the orthogonal
%   polynomials for the set of points (x,y). P is the orthogonal
%   matrix which contains as columns the values p_i(x) of the
%   orthogonal polynomials, and b contains the coefficients for the
%   expansion P(x)=b_1 p_1(x) + ... + b_m p_m(x).

m=length(x); x=x(:); y=y(:);
alpha(1)=sum(x)/m;
p1=ones(size(x)); p2=x-alpha(1);
P=[p1,p2];
b(1)=p1'*y/norm(p1)^2; b(2)=p2'*y/norm(p2)^2;
for k=1:m-2
  p0=p1; p1=p2;
  alpha(k+1)=p1'*(x.*p1)/norm(p1)^2;
  beta(k)=(norm(p1)/norm(p0))^2;
  p2=(x-alpha(k+1)).*p1-beta(k)*p0;
  P=[P,p2];
  b(k+2)=p2'*y/norm(p2)^2;
end
b=b(:);