工作中需要用到polylog方程。但是在matlab里面没有这个方程的定义。找了很久，只发现有polylog的二阶定义

dilog。相比之下，mathematics里面有polylog的函数，应用起来很方便。因为自己用matlab比较顺手，所以还是希望能有polylog函数在matlab里面的定义。后来在另一个同事中找到一个matlab的近似定义，我比较了2阶情况。最大误差在1e-3。或许对其他人也有帮助。

%This function will calculate the polylog g_s(z) for -inf<z<1,

%given that z is a real number.  The cutoff for switching to the

%asymptotic expansion limits the precision to 1e-4 in this regime.

%The desired precision for -1<z<1 can be specified with "error".

function out = polylog(z,s,error)

if s>1

   %perform piece-wise calculation of polylog g_s(z) with spec'd error

   j = 1;              %vector index

   maxj = length(z);   %input vector length

   g = zeros(size(z)); %initialize output

   while j<=maxj

   if abs(z(j))<=1        %series definition

   zn = 1;             %z^n term

   epsilon = 1;        %error

   n = 1;              %sum index

   while abs(epsilon)>error

       zn = zn.*z(j);

       gn = zn./n.^s;      %correction for current step

       g(j) = g(j) + gn;

       epsilon = gn./g(j);

       n = n+1;

   end

   elseif z(j)<-500             %asymptotic expansion, Eq. 15, [Pathria, R.K., Statistical Mechanics 2nd Ed., 510 (1972)]

       xsi = log(-z(j));

       g1 = xsi.^s./gamma(s+1);

       g(j) = g1.*(1 + s*(s-1)*pi^2/6./xsi.^2);

       g(j) = g(j) + g1.*(s*(s-1)*(s-2)*(s-3)*7*pi^4/360./xsi.^4);

       g(j) = g(j) + g1.*(s*(s-1)*(s-2)*(s-3)*(s-4)*(s-5)*31*pi^6/15120./xsi.^6);

       %epsilon = g1.*(s*(s-1)*(s-2)*(s-3)*(s-4)*(s-5)*31*pi^6/15120./xsi.^6)

       g(j) = -g(j);

   else                       %numerical integration, Eq. 2, pg. 508, Pathria

       g(j) = 1/gamma(s);

       x = linspace(0,30,1000);

       f = x.^(s-1)./(exp(x)./(-z(j))+1);

       g(j) = -g(j).*trapz(x,f);

   end

   j=j+1;

   end

   %output

   out = g;

   %References

   %asymptotic expansion, Eq. 15, [Pathria, R.K., Statistical Mechanics 2nd

   %Ed., 510 (1972)]

   %f_s(z) = xsi^s/gamma(s+1)*(1 + 2*s*sum_{j=1,3,5...}(s-1)...(s-j)*(1-1/2^j)*zeta(j+1)/xsi^(j+1))

else

   out=-(polylog(z,s+1,error)-polylog(z+z*0.001,s+1,error))/0.001;

end;