I am trying to find the distribution function of a random variable $Z = X-Y$ with $X$ following a $\text{Hypo}(\mu_1,...,\mu_n)$ and $Y$ following a $\text{Hypo}(\lambda_1,...,\lambda_m)$. Then $X$ and $Y$ both have the support region $[0,\infty]$. Additionally, I assumed that all parameters are distinct such that we can use the PDFs:
\begin{equation} f_X(x) = \sum_{k=1}^{n}l_{k}(0)\mu_ke^{-\mu_kx}, \; f_Y(y) = \sum_{k=1}^{m}l_{k}(0)\lambda_ke^{-\lambda_ky} \end{equation}
with $l_1(x),...,l_n(x)$ are the Lagrange basis polynomials.
I tried to use the convolution theorem for probability theorem, such that:
\begin{equation} f_Z(z) = \int_{0}^{\infty}f_{X}(z+y)f_Y(y)\text{d}y \end{equation}
but I end up obtaining something of the form:
\begin{equation} \sum_{k_1 = 1}^{n}\sum_{k_2=1}^{m}\frac{l_{k_1}(0)l_{k_2}(0)}{\mu_{k_1} + \lambda_{k_2}} \mu_{k_1} \lambda_{k_2} e^{-\mu_{k_1}z} \end{equation}
which is a divergent when you take the integral over the support region of Z, which I believe is $[-\infty, \infty]$.
Is there somebody more familiar with this topic and willing to help me?
In my opinion you computed the convolution in the wrong way, since you didn't account for z positive or negative. I would do it like this:
$f_Z(z) = \int_{\mathbb{R}}f_X (x) f_Y(x-z) dx $
When $z <0$ the integral is evaluated in $(0,\infty)$, while for $z \ge 0$ in $(z,\infty)$. In the end I get the following density for $Z$:
$$ f_Z(z) = \sum_{k=1}^{n}\sum_{s=1}^{m} \frac{l_s(0)l_k(0)\lambda_s\mu_k}{\lambda_s+\mu_k} (\mathbb{I}_{\{ z <0 \}} + e^{-(\lambda_s+\mu_k)z}\mathbb{I}_{\{ z \ge 0 \}})$$