function [P,u,v] = forsythe(xi,q,m)
%
%  Name: forsythe
%
%  Usage: [P,u,v] = forsythe(xi,q,m)
%
%  Description: Generates orthogonal Forsythe
%		polynomials over the domain 'xi'
%		with the weights 'q' of order 'm'
%
%		Returns matrix 'P', which has 
%		size (L x m+1). (L = length(xi))
%		and normalization factors 'u' and
%		'v'
%

%  Version SWD960507
%
%  Author(s): Scott W. Doebling, doebling@lanl.gov
%
%  Institution: Los Alamos National Laboratory,
%		ESA-EA
%
%  Copyright (c) 1996, The Regents of the University
%     	of California. (See top level of code for
%	full copyright statement.)

L = length(xi);
P = zeros(L,m+1);
u = zeros(m+1,1);
v = zeros(m+1,1);
      
%
% Define zero-th order poly.
%

P(:,1) = ones(L,1);

%
% Define first order poly.
%
      
d(1) = sum( P(:,1).^2 .* q.^2);
u(2) = sum( xi .* P(:,1).^2 .* q.^2) / d(1);

P(:,2) = (xi - u(2)) .* P(:,1);

%
% Define higher-order polys.
%

for k = 3:(m+1),

   d(k-1) = sum( P(:,k-1).^2 .* q.^2);

   v(k-1) = sum( xi .* P(:,k-1) .* P(:,k-2) .* q.^2) / d(k-2);

   u(k) = sum( xi .* P(:,k-1).^2 .* q.^2) / d(k-1);
   
   P(:,k) = (xi - u(k)) .* P(:,k-1) - v(k-1) .* P(:,k-2);

end
      
d(m+1) = sum( P(:,m+1).^2 .* q.^2);

return