Hi,
Indeed, in Mathematica version 12.2.0 I managed to get the Inverse Laplace transform for the k-th coefficient using Apart[], as suggested by larawag. I chose k to be from 1 to 40, albeit slow it eventually gives the correct result.
However this trick is not enough for solving the problem if you use older versions of Mathematica. I suggest following workaround.
I found this useful InverseLaplaceTransform from the book "Tables of Integral Transforms: Based, in Part, on Notes Left by Harry Bateman":
Using it I wrote a small code:
In[1]:= $Version
Out[1]= "11.3.0 for Mac OS X x86 (64-bit) (March 7, 2018)"
In[2]:= ILT[f_, s_, t_] /; Head[f] === Plus := ILT[#, s, t] & /@ f;
ILT[a_ f_, s_, t_] /; FreeQ[a, s] := a*ILT[f, s, t];
ILT[s_^mu_ *(s_ + b_)^nu_, s_, t_] := Simplify[
ConditionalExpression[
1/Gamma[-nu - mu] t^(-nu - mu - 1)
Hypergeometric1F1[-nu, -nu - mu, -b t], nu + mu < 0]];
ILT[f_, s_, t_] := InverseLaplaceTransform[f, s, t];
In[6]:= nmax = 7;
In[7]:= F[
s_] := (((1 + s)^2 + 3*s^(3/2)*(1 + s)^(1/2)*th - 4*s*(1 + s)*th +
3*s^(1/2)*(1 + s)^(3/2)*th +
s^2*th^2))/(8 (s^(3/2)*(1 + s)^(1/2)*((1 + s)^(1/2) +
s^(1/2)*th)^3));
In[8]:= ser0 = Simplify[Series[F[s], {th, 0, nmax}]];
In[9]:= coeffs0 = Table[SeriesCoefficient[ser0, n], {n, 0, nmax}];
In[10]:= coeffs1 = ILT[Apart[#], s, t] & /@ coeffs0;
In[11]:= Simplify[coeffs1]
Out[11]= {Sqrt[t]/(4 Sqrt[\[Pi]]), (
E^-t (3 Sqrt[t] - 2 Sqrt[\[Pi]] Erfi[Sqrt[t]]))/(4 Sqrt[\[Pi]]), (
E^-t (-8 (-3 + E^t) Sqrt[t] +
Sqrt[\[Pi]] (-7 + 8 t) Erfi[Sqrt[t]]))/(8 Sqrt[\[Pi]]), (
E^-t ((23 - 12 E^t - 10 t) Sqrt[t] +
6 Sqrt[\[Pi]] (-1 + 2 t) Erfi[Sqrt[t]]))/(4 Sqrt[\[Pi]]), (
E^-t (2 Sqrt[t] (80 (3 - 2 t) + 9 E^t (-15 + 4 t)) -
9 Sqrt[\[Pi]] (11 - 34 t + 8 t^2) Erfi[Sqrt[t]]))/(
48 Sqrt[\[Pi]]), (
E^-t (2 Sqrt[t] (354 - 376 t + 56 t^2 + 45 E^t (-5 + 2 t)) -
45 Sqrt[\[Pi]] (3 - 12 t + 4 t^2) Erfi[Sqrt[t]]))/(
48 Sqrt[\[Pi]]), (1/(320 Sqrt[\[Pi]]))
E^-t (-2 Sqrt[
t] (-224 (15 - 20 t + 4 t^2) + 5 E^t (439 - 294 t + 32 t^2)) +
5 Sqrt[\[Pi]] (-225 + 1140 t - 620 t^2 + 64 t^3) Erfi[Sqrt[
t]]), (1/(120 Sqrt[\[Pi]]))
E^-t (-2 Sqrt[
t] (35 E^t (33 - 28 t + 4 t^2) +
3 (-555 + 950 t - 316 t^2 + 24 t^3)) +
35 Sqrt[\[Pi]] (-15 + 90 t - 60 t^2 + 8 t^3) Erfi[Sqrt[t]])}
The first two lines in In[2] are just to ensure linearity for the function ILT, which will be used to find the Inverse Laplace transform.
I also put nmax=40 and on my machine it took only couple of seconds to get the final result.
Hope this was useful,
Hrach