function [yy,pp]=OrthogonalInterpolation(alpha,beta,b,xx)
% ORTHOGONALINTERPOLATION interpolation using orthogonal polynomials
%   [yy,pp]=OrthogonalInterpolation(alpha,beta,b,xx) evaluates the
%   precomputed orthogonal polynomials defined by alpha, beta, b from
%   OrthogonalPolynomialCoefficients(x,y) for the values xx and
%   returns their values in pp and the approximations yy(:,i)=b_1
%   p_1(xx) + ... + b_i p_i(xx)

m=length(b);
p1=ones(size(xx)); p2=xx-alpha(1);
pp=[p1,p2];  
yy=[p1*b(1) p1*b(1)+p2*b(2)];
for i=2:m-1
  p0=p1; p1=p2;
  p2=(xx-alpha(i)).*p1-beta(i-1)*p0;
  pp=[pp,p2];
  yy=[yy yy(:,i)+p2*b(i+1)];
end