We develop an algorithm for computing the (exact) cumulative distribution function of the time-dependent workload in a piecewise-stationary Mt/G(t)/1 queue with a work-conserving service discipline and general service-time distributions, where service times are determined at arrival instants. The t subscripts indicate that the arrival rate and the general service-time distribution may change with time, but we allow changes only at finitely many time points. The algorithm is based on numerical transform inversion, using the classical Takacs double-transform of the transient workload in a M/G/1 queue recursively over the successive stationary intervals. In particular, we apply our recently developed Fourier-series-based inversion algorithms for two-dimensional transforms and nested one-dimensional transforms. We also do additional work to greatly speed up the computation while tightly controlling the error. As a consequence, the computation time grows only quadratically with the number of intervals. The algorithm is effective for ten or fewer intervals, where the intervals may have unlimited and possibly unequal lengths, typically running in at most a few minutes and maintaining high accuracy. We have also demonstrated that the algorithm can solve a 21-interval example with 7-to-10-digit accuracy in about half an hour. Models with only a few intervals are useful to study overload control strategies.
