How to quickly calculate the following function？

Calculating the sum seems to be more stable. For cases, in which the integral method is successful, the sum is 10 to 200 times faster for the given test data:

function test_Lerch
% See: https://people.math.sc.edu/Burkardt/py_src/polpak/lerch_values.py
% and: https://people.sc.fsu.edu/~jburkardt/m_src/test_values/test_values.html
a_vec = [0.0E+00, ...
   0.0E+00, ...
   0.0E+00, ...
   1.0E+00, ...
   1.0E+00, ...
   1.0E+00, ...
   2.0E+00, ...
   2.0E+00, ...
   2.0E+00, ...
   3.0E+00, ...
   3.0E+00, ...
   3.0E+00];
s_vec = [2, 3, 10, ...
   2, 3, 10, ...
   2, 3, 10, ...
   2, 3, 10];
z_vec = [0.1000000000000000E+01, ...
   0.1000000000000000E+01, ...
   0.1000000000000000E+01, ...
   0.5000000000000000E+00, ...
   0.5000000000000000E+00, ...
   0.5000000000000000E+00, ...
   0.3333333333333333E+00, ...
   0.3333333333333333E+00, ...
   0.3333333333333333E+00, ...
   0.1000000000000000E+00, ...
   0.1000000000000000E+00, ...
   0.1000000000000000E+00];
f_vec = [ ...
   0.1644934066848226E+01, ...
   0.1202056903159594E+01, ...
   0.1000994575127818E+01, ...
   0.1164481052930025E+01, ...
   0.1074426387216080E+01, ...
   0.1000492641212014E+01, ...
   0.2959190697935714E+00, ...
   0.1394507503935608E+00, ...
   0.9823175058446061E-03, ...
   0.1177910993911311E+00, ...
   0.3868447922298962E-01, ...
   0.1703149614186634E-04];
for k = 1:numel(z_vec)
   z = z_vec(k);
   s = s_vec(k);
   a = a_vec(k);
   f = f_vec(k);
   
   fprintf('\nk: %d\n', k)
   tic
   for k = 1:1
      y1 = Lerch_fun_integral(z,s,a);
   end
   toc
   
   tic
   for k = 1:1
      y2 = Lerch_fun_sum(z,s,a);
   end
   toc
   
   fprintf('Ref: %.16g\n', f);
   fprintf('Int: %.16g  delta: %16g\n', y1, abs(y1 - f));
   fprintf('Sum: %.16g  delta: %16g\n', y2, abs(y2 - f));
end
end
function out = Lerch_fun_integral(z,s,a)
f=@(t)((t.^(s-1)).*exp(-a*t))./(1-z.*exp(-t));
out=(1/gamma(s))*(integral(f,0,Inf));
end
function out = Lerch_fun_sum(z, s, a)
out = 0;
if z <= 0
   return
end
lim = 1e-16;
k   = 0;
z_k = 1;
term = Inf;
while abs(term) > lim * (1.0 + abs(out))
   if a + k ~= 0
      term = z_k / (a + k)^s;
      out  = out + term;
   end
   k   = k + 1;
   z_k = z_k * z;
end
end