After taking the inverse Laplace transform of a function, I wrote a Maple command to evaluate the result at some specified values of t (for example at t=0), it returned solution in equation (2) . I expected single numerical value for each value of t, but , it gave the result in form of linear combinations of delta(2,0), delta(3,0), delta(4,0), etc. See below:

> restart;

> with(Physics);

> Setup(mathematicalnotation = true); complex*number;

[mathematicalnotation = true]

complex number

> k[1] := 139.39; c[1] := .2725; mu[1] := 5.9410; gamm1 := 2;

> P[0] := 1; M := 70; g := 9.8; v := 6.5; L := 6;

> P1 := P[0]*M*g*(-L*s/v)^k/factorial(k);

> Y[1](0), Y[1](1), Y[1](2), Y[1](3) := 0., a, 0., c;

> Y[2](0), Y[2](1), Y[2](2), Y[2](3) := 0., b, 0., d;

> B[k-r] := KroneckerDelta[k, r]-(3/2)*KroneckerDelta[k, r+1]+(3/4)*KroneckerDelta[k, r+2]-(1/8)*KroneckerDelta[k, r+3];

> B1 := -(3/2)*KroneckerDelta[k, r]+(3/2)*KroneckerDelta[k, r+1]-(3/8)*KroneckerDelta[k, r+2];

> B2 := (3/2)*KroneckerDelta[k, r]-(3/4)*KroneckerDelta[k, r+1];

> Q[k-r] := KroneckerDelta[k, r]-(1/2)*KroneckerDelta[k, r+1];

> A1 := sum(B[k-r]*(r+1)*(r+2)*(r+3)*(r+4)*Y[1](r+4), r = 0 .. k)+2*(sum(B1*(r+1)*(r+2)*(r+3)*Y[1](r+3), r = 0 .. k))+sum(B2*(r+1)*(r+2)*Y[1](r+2), r = 0 .. k)+k[1]*(Y[1](k)-Y[2](k))+c[1]*s*(Y[1](k)-Y[2](k))-mu[1]*s*(k+1)*(k+2)*(Y[1](k+2)-Y[2](k+2))+sum(gamm1*s^2*Q[k-r]*Y[1](r), r = 0 .. k) = P1;

> A2 := sum(B[k-r]*(r+1)*(r+2)*(r+3)*(r+4)*Y[2](r+4), r = 0 .. k)+2*(sum(B1*(r+1)*(r+2)*(r+3)*Y[2](r+3), r = 0 .. k))+sum(B2*(r+1)*(r+2)*Y[2](r+2), r = 0 .. k)+k[1]*(Y[2](k)-Y[1](k))+c[1]*s*(Y[2](k)-Y[1](k))-mu[1]*s*(k+1)*(k+2)*(Y[2](k+2)-Y[1](k+2))+sum(gamm1*s^2*Q[k-r]*Y[2](r), r = 0 .. k) = 0;

> m := 12;

> for i from 0 to m do Y[1](i+4) := collect(solve(eval(A1, k = i), Y[1](i+4)), {a, b, c, d}); Y[2](i+4) := collect(solve(eval(A2, k = i), Y[2](i+4)), {a, b, c, d}) end do;

>

> A3 := sum(Y[1](k), k = 0 .. 8) = 0;

> A4 := sum(k*(k-1)*Y[1](k), k = 0 .. 8) = 0;

> A5 := sum(Y[2](k), k = 0 .. 8) = 0;

> A6 := sum(k*(k-1)*Y[2](k), k = 0 .. 8) = 0;

> sol := solve({A3, A4, A5, A6}, {a, b, c, d});

> Y[1] := eval(sum(Y[1](k)*xi^k, k = 0 .. 8), sol);

>

> with(inttrans);

> y[1] := invlaplace(Y[1], s, t);

> y[2] := subs(I = 0, y[1]);

>

> for xi from 0 by .1 to 1 do y[xi] := eval(eval(y[2], {t = 0})) end do;

y[0] := 0.

y[0.1] := -335.8210943 - 0.2853134161 delta(1, 0)

- 0.003189135039 delta(2, 0) + 3.215933203 delta(0)

-8 -9

- 2.101888941 10 delta(3, 0) + 8.278917357 10 delta(4, 0)

y[0.2] := -616.5812812 - 0.7286587659 delta(1, 0)

+ 0.003338719654 delta(2, 0) + 2.597336701 delta(0)

-7

- 0.000001136393191 delta(3, 0) + 2.587772223 10 delta(4, 0)

y[0.3] := -787.6445525 - 1.474810997 delta(1, 0)

+ 0.02853884261 delta(2, 0) - 5.517159103 delta(0)

- 0.00001031142443 delta(3, 0) + 0.000001812552944 delta(4, 0)

y[0.4] := -811.1848471 - 2.589947077 delta(1, 0)

+ 0.07684191276 delta(2, 0) - 23.27219517 delta(0)

- 0.00004426296248 delta(3, 0) + 0.000006660566248 delta(4, 0)

y[0.5] := -682.8521581 - 3.984266525 delta(1, 0)

+ 0.1440502970 delta(2, 0) - 49.60730104 delta(0)

- 0.0001229043255 delta(3, 0) + 0.00001666149585 delta(4, 0)

y[0.6] := -441.5562325 - 5.350681306 delta(1, 0)

+ 0.2144868917 delta(2, 0) - 78.78970103 delta(0)

- 0.0002502887215 delta(3, 0) + 0.00003148839539 delta(4, 0)

y[0.7] := -167.0511892 - 6.158350095 delta(1, 0)

+ 0.2612340304 delta(2, 0) - 100.0143657 delta(0)

- 0.0003906180172 delta(3, 0) + 0.00004649544839 delta(4, 0)

y[0.8] := 39.3465739 - 5.754055618 delta(1, 0)

+ 0.2517612354 delta(2, 0) - 99.27848038 delta(0)

- 0.0004565738656 delta(3, 0) + 0.00005214660880 delta(4, 0)

y[0.9] := 94.9253902 - 3.638647825 delta(1, 0)

+ 0.1617018641 delta(2, 0) - 65.09861522 delta(0)

- 0.0003369182899 delta(3, 0) + 0.00003737844304 delta(4, 0)

-7

y[1.0] := -0.00034023 + 5.342230828 10 delta(1, 0)

-8

- 2.395810153 10 delta(2, 0) + 0.000005214637385 delta(0)

-9 -11

- 1.077399852 10 delta(3, 0) + 1.681923081 10 delta(4, 0)