The Laplace transform is an invaluable tool for analyzing Markov processes. In continuous-time Markov chains, the distribution at time $t$, expressed as $\pi(t) = \pi_0 \exp(tQ)$, is more complex than its Laplace transform $\hat{\pi}(s) = \int_0^\infty \pi_0 \exp(tQ)e^{-st} \,dt = \pi_0(Q - sI)^{-1}$, which requires only a matrix inverse rather than a matrix exponential \[5]. This analogy extends to more complex models such as Markov-modulated fluid models $(X, \phi)$, where the first-return matrix $[\Psi(t)]_{ij} = P[\phi(\tau) = j, \tau < t \mid X(0) = 0, \phi(0) = i]$, with $\tau = \min(t > 0 \mid X(t) = 0)$, is hard to compute directly; however, efficient numerical algorithms exist for computing its Laplace transform \[2]. While working with Laplace transforms is convenient, inverting them to compute $f(t)$ from $\hat{f}(s)$ is numerically challenging. Unlike Fourier transforms, inverse Laplace transforms are ill-posed—small perturbations in $\hat{f}(s)$ can lead to large changes in $f(t)$. Many inversion methods have been developed, including the Abate–Whitt framework \[1], which approximates $f(t) \approx \sum_{i=1}^n \frac{w_i}{t} \hat{f}\left(\frac{\beta_i}{t}\right)$ using specified weights $w_i$ and nodes $\beta_i \in \mathbb{C}$, encompassing classical techniques such as Euler’s method. More recently, the CME method \[3] was introduced, which approximates the Dirac delta $\delta_1(s)$ with a weighted sum of exponentials. Many of these approaches are designed to handle irregular functions with discontinuities, but matrix-exponential distributions, which often arise in applications, are much more regular—suggesting a more favorable setting for inversion. This talk focuses on the specific task of computing inverse Laplace transforms for such probability-related functions. We demonstrate that, in this context, the inversion error can be bounded by the approximation error $\max_{y \in Y} \left| e^y - \sum_{i=1}^n \frac{w_i}{\beta_i - y} \right|$, over a region $Y \subset \mathbb{C}$ that includes the eigenvalues of $tQ$ (in the Markov chain case) or an analogous quantity for fluid models. A convenient proxy for this region can be derived using the optimal uniformization rate $\lambda = \max_i |Q_{ii}|$, which ensures the eigenvalues of $Q$ lie within a circle centered at $-\lambda$ with radius $\lambda$. We revisit a lesser-known method within the Abate–Whitt framework, due to Zakian, which approximates the exponential function using a Padé approximant. We then highlight a more recent numerical method—the AAA (adaptive Antoulas–Anderson) approximation algorithm \[4]—which yields high-quality rational approximants and better parameter families for use within the Abate–Whitt framework. Our analysis further reveals that the expected accuracy of these methods depends on both $\lambda$ and the evaluation point $t$; problems requiring approximation of the exponential over a broader region tend to be more challenging and demand a greater number of nodes.