function S=sonine(p,l)
%
% Sonine Polynomial 
%    S = sonine(p,l) finds the coefficients of the sonine polynomial S(p,l)
%    of degree p. S is a row vector of length p+1 containing the polynomial
%    coefficients in descending powers, S(1)*X^p + S(2)*X^(p-1) +...+
%    S(p)*X + S(p+1). 
%
%    Use sonine(p,d/2-1) for the d dimensional expansion.

n=p:-1:0;

S=gamma(p+l+1)*(-1).^n./(gamma(n+l+1).*factorial(p-n).*factorial(n));