Hypoexponential distribution is stuck because of $\lambda_i \neq \lambda_j$? The definition of the hypoexponential distribution (HD) requires that:
$$f(x)=\sum_i^d \left(\prod_{j=1,i\neq j}^{d}\frac{\lambda_j}{\lambda_j-\lambda_i}\right)\lambda_i e^{(-\lambda_ix)},\quad x>0
$$
is well-defined, which imposes that all $\lambda_i \neq \lambda_j$. As such, it seems to impose a significant limitation on the underlying exponentially distributed $X_i$'s in that if, say $X_i$ and $X_j$, have $\lambda_i=\lambda_j$, the distribution is no longer HD, neither can be Erlang, which requires that all $\lambda_i$ are the same!
I'm sure I'm missing something!
A second question: I tried to plot the HD using the syntax below for $\lambda = .1:20$ but not sure if it's the right syntax. Maybe you have experience this this R package ir can suggest another one for working with HD (generalized Erlang)?
library(sdprisk)
x <- seq(0.1, 20, 0.05)
y <- dhypoexp(x, rate = r, log = FALSE)
plot(x, y)




The definition of the hypoexponential distribution (HD) requires that:
$$f(x)=\sum_i^d \left(\prod_{j=1,i\neq j}^{d}\frac{\lambda_j}{\lambda_j-\lambda_i}\right)\lambda_i e^{(-\lambda_ix)},\quad x>0
$$

that expression is only a special case when $\lambda_i \neq \lambda_j$. It is not a requirement for the hypoexponential distribution in general.
The more general hypo-exponential distribution can be expressed as a phase-type distribution
$$f(x) = -\boldsymbol{\alpha}e^{x \boldsymbol{\Theta}} \boldsymbol{\Theta} \mathbf{1} $$
With $$\boldsymbol{\alpha} = (1,0,0,\dots,0,0)$$ and
$$\boldsymbol{\Theta} = \begin{bmatrix} 
-\lambda_1 & \lambda_1 & 0 & \dots & 0 & 0 \\
0 & -\lambda_2 &  \lambda_2 & \dots & 0 & 0 \\
\vdots & \ddots & \ddots & \ddots & \ddots & \vdots\\
0 & 0 & \ddots & -\lambda_{d-2} & \lambda_{d-2} & 0\\
0 & 0 & \dots & 0 & -\lambda_{d-1} & \lambda_{d-1}\\
0 & 0 & \dots & 0 & 0 & -\lambda_{d}\\
\end{bmatrix}$$
That involves the exponentiation of a matrix. And that can be approximated with a sum using a Taylor series.
You can see the matrix as modelling a sort of Markov chain process with non-discrete time steps. A sum of exponential distributed variables is like waiting for several consecutive transitions whose waiting time is each exponential distributed. Those transitions relate to the Markov chain.
Not only when $\lambda_i \neq \lambda_j$ but also when $\lambda_i = \lambda_j$ then the formula can give problems as demonstrated in this question: Why my cdf of the convolution of n exponential distribution is not in the range(0,1)?
