function JacPols = JacobiPolynomials(alpha,beta,max_deg)

% computes normalized Jacobi polynomials of type (alpha,beta) upto the maximum degree = max_deg, the are return in a cell array

a_n = @(n) (2 / (2 * n + alpha + beta)) * sqrt((n * (n + alpha + beta) * (n + alpha) * (n + beta)) / ((2 * n + alpha + beta - 1) * (2 * n + alpha + beta + 1)));
b_n = @(n) -(alpha ^ 2 - beta ^ 2) / ((2 * n + alpha + beta) * (2 * n + alpha + beta + 2));

J0 = sqrt((2 ^ (-alpha - beta - 1)) * ((gamma(alpha + beta + 2)) / (gamma(alpha + 1) * gamma(beta + 1))));
J1 = 0.5 * J0 * sqrt((alpha + beta + 3) / ((alpha + 1) * (beta + 1))) * [(alpha + beta + 2) (alpha - beta)];
JacPols = cell(max_deg + 1,1);
JacPols{1} = J0;
if max_deg > 0
	JacPols{2} = J1;
end
if max_deg > 1
	for n = 1:(max_deg - 1)
		JacPols{n + 2} = ([JacPols{n + 1} 0] - b_n(n) * [0 JacPols{n + 1}] - a_n(n) * [0 0 JacPols{n}]) / a_n(n + 1);
	end
end

end
