matrix $\mH$ and any constant $c>0$, we denote by $$\begin{aligned}
(\Lambda_{c}^{+}(\mH), \Lambda_{c}^{-}(\mH),\mG_{c}^{+}(\mH), \mG_{c}^{-}(\mH))\end{aligned}$$ the unique quadruple constituting the classical Wiener-Hopf factorization (cf. [@Barlow1980]) corresponding to $\mH$ with killing rate $c$. In order to proceed, we let $c_{k}=q_{k}+c,\quad k\geq 1$.

We are now ready to describe the algorithm to compute $q_{1}^{-1}\cdots q_{n}^{-1}\widehat{\Pi}_{c}^{+}(i,j;q_1,\ldots,q_n)$.

1.  **Compute the first diagonal:** for $k=1,\ldots,n+1$, compute $$\begin{aligned}
    {\widetilde}{\Lambda}_{c,k-1,k-1}^{+}=\Lambda_{c_{k}}^{+}(\mG_{k}),\end{aligned}$$ using the diagonalization method in [@rogers_shi_1994].

2.  **Compute the second diagonal:** for $k=1,\ldots,n$, solve the following linear system $$\begin{aligned}
    q_{k}\mI^{+}+\mB_{k}{\widetilde}{\Lambda}_{c,k-1,k}^{+}&=\mV^{+}{\widetilde}{\mG}_{c,k-1,k}^{+},\\
    [\mD_{k}-c_{k}\mI^{-}]{\widetilde}{\Lambda}_{c,k-1,k}^{+}+q_{k}{\widetilde}{\Lambda}_{c,kk}^{+}&=\mV^{-}{\widetilde}{\Lambda}_{c,k-1,k-1}^{+}{\widetilde}{\mG}_{c,k-1,k}^{+}+\mV^{-}{\widetilde}{\Lambda}_{c,k-1,k}^{+}{\widetilde}{\mG}_{c,kk}^{+},\end{aligned}$$ for ${\widetilde}{\Lambda}_{c,k-1,k}^{+}$ and ${\widetilde}{\mG}_{c,k-1,k}^{+}$.

3.  **Compute the other diagonals:** for $r=2,\ldots,n$, $k=0,\ldots,n-r$, solve the linear system $$\begin{aligned}
    \mB_{k+1}{\widetilde}{\Lambda}_{c,k,k+r}^{+}&=\mV^{+}{\widetilde}{\mG}_{c,k,k+r}^{+},\\
    [\mD_{k+1}-c_{k+1}\mI^{-}]{\widetilde}{\Lambda}_{c,k,k+r}^{+}+q_{k+1}{\widetilde}{\Lambda}_{c,k+1,k+r}^{+}&=\mV^{-}\sum_{j=0}^{r}{\widetilde}{\Lambda}_{c,k,k+j}^{+}{\widetilde}{\mG}_{c,k+j,k+r}^{+},\end{aligned}$$ for ${\widetilde}{\Lambda}_{c,k,k+r}^{+}$ and ${\widetilde}{\mG}_{c,k,k+r}^{+}$.

4.  **Compute** $$\begin{aligned}
    P^{+}(q_{1},\ldots,q_{n}):=q_{1}^{-1}\cdots q_{n}^{-1}\widehat{\Pi}_{c}^{+}(i,j;q_1,\ldots,q_n)=q_{1}^{-1}\cdots q_{n}^{-1}\sum_{\ell=0}^{n}\,{\widetilde}{\Lambda}_{c,0\ell}.\end{aligned}$$ for $q_1,\ldots,q_n>0$.

5.  **Compute the approximate inverse Laplace transform of $P^{+}(q_{1},\ldots,q_{n})$:** use the method discussed in Section \[subsubsec:specialcase\].

If $|\mathbf{E}^{+}|=|\mathbf{E}^{-}|=1$, then the matrices in Steps 1-3 become numbers. Step 1 reduces to solving $n+1$ quadratic equations for a root in $[0,1]$. In Step 2 and 3, for each loop, the system reduces to a system of two linear equations of two unknowns. Moreover, in this case, $P^{+}$ has a closed-form representation for $q_1,\ldots,q_n>0$, and hence, for any $q_1,\ldots,q_n\in\bC^+$, as mentioned in the previous section. This allows to use general 