function Ps = JacobiPolRek(n,alpha,beta,x)
% Rekursive Berechnung der Jacobi-Polynome
% P_n = P(n+1)  (Struktur)
% Aufruf  P_n = JacobiPolRek(n,alpha,beta,x)
% typische Werte und Bedeutung der Parameter vgl. Buch
% Bsp.:
% phi = linspace(0,pi);
% P3=JacobiPolRek(3,2,1/2,cos(phi));
% P3
% 
% P3 = 
% 
%         P: [4x100 double]
%         n: 3
%     alpha: 2
%      beta: 0.5000
%         x: [1x100 double]

% Vers. 1.2 14.08.2016
% Copyright Wolfgang Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016


n=n+1;
P = ones(size(x));
P = repmat(P,n,1);
%P(1) = 1;.x
P(2,:) = 1/2*(alpha-beta+(alpha+beta+2).*x);
for k=3:n
    a(1) = 2*k*(k+alpha+beta)*(2*k-2+alpha+beta);
    a(2) = (2*k-1+alpha+beta)*(alpha^2-beta^2);
    a(3) = (2*k-2+alpha+beta)*(2*k-1+alpha+beta).*(2*k+alpha+beta);
    a(4) = 2*(k-1+alpha)*(k-1+beta)*(2*k+alpha+beta);
    %
    P(k,:) = ((a(2)+a(3) .* x).*P(k-1,:) - a(4).*P(k-2,:))./a(1);
end
Ps.P = P;
Ps.n = n-1;
Ps.alpha = alpha;
Ps.beta = beta;
Ps.x = x;
