Problem with Inverse Laplace Transform in Mathematica

Question

I am trying to calculate the following inverse Laplace transform

$$ f(t) = \mathcal{L}^{-1}\left[ \left(\frac{1}{s+a}\right)^{n+1}\left( \frac{1}{s+b}\right)^{n} \right](t) $$

with $s$ being the variable, $a,b >0$ and $n\in N$. I know I can use the properties of the inverse transform to calculate this from the convolution

$$ \mathcal{L}^{-1}\left[ \left(\frac{1}{s+a}\right)^{n+1}\right]*\mathcal{L}^{-1}\left[\left( \frac{1}{s+b}\right)^{n} \right](t) $$

and the fact that:

$$ \mathcal{L}^{-1}\left[ \left(\frac{1}{s+a}\right)^{n+1}\right] = \frac{t^n e^{-at}}{n!}\Theta(t) $$

Combining this together this gives for the final result of the inverse transform:

$$ f(t) = \frac{\sqrt{\pi t (b-a)}}{2 \Gamma(n+1)} e^{-\frac{1}{2}(a+b)t} \left[I_{n-1/2}\left( \frac{1}{2}(b - a) t\right)+I_{n+1/2}\left( \frac{1}{2} (b - a) t \right) \right] , $$

where $I_{\nu}(x)$ is the Modified Bessel function of the first kind. However, if I try to check whether or not this is correct by taking the Inverse transform I get an expression that only gives the initial form if $b>a$, but not otherwise

I cannot understand what the problem is, or if I am missing something important.

Here is the Mathematica code I am using:

fun1 = 1/(s + a)^(NN + 1)

fun2 = 1/(s + b)^NN 

fun = fun1*fun2

baseTransform = 1/(s + x)^M

baseInverse = InverseLaplaceTransform[baseTransform, s, t1]

fN = baseInverse /. {M -> (NN + 1), x -> a, t1 -> t - tau}

gN = baseInverse /. {M -> NN, x -> b, t1 -> tau}

sol = Refine[Integrate[fN*gN, {tau, 0, t}], {NN \[Epsilon] Integer, NN > 0, t \[Epsilon] Real, t > 0}] 

funCheck1 = Simplify[LaplaceTransform[sol, t, s],Assumptions -> {a < b}]

funCheck2 = Simplify[LaplaceTransform[sol, t, s],Assumptions -> {a > b}]

Answer 1

If we use:

L1 = InverseLaplaceTransform[1/(s + a)^(n + 1), s, t]; 
L2 = InverseLaplaceTransform[1/(s + b)^n, s, t];
L12 = FullSimplify[Integrate[(L1 /. t -> (t - \[Tau]))*(L2 /. t -> \[Tau]), {\[Tau], 0,t}, Assumptions -> n > 0], t > 0]

f[a_, b_, n_, t_] := InverseLaplaceTransform[1/(s + a)^(n + 1) 1/(s + b)^n, s, t] // N; 
g[a_, b_, n_, t_] := ((-a + b)^(1/2 - n) E^(-(1/2) (a + b) t) Sqrt[\[Pi]] t^(
1/2 + n) (BesselI[-(1/2) + n, 1/2 (-a + b) t] + BesselI[1/2 + n, 1/2 (-a + b) t]))/(2 Gamma[1 + n]) // N

Then we can check and works fine for a > b and b < a:

(Table[g[a, b, n, 1] - f[a, b, n, 1] /. Indeterminate -> 0, {a, -5, 5}, {b, -5, 5}, {n, -5, 5}] // Quiet // Chop) // MatrixForm(* for t = 1 *)

(*ALL Zeros *)