Amacroscopic two-fluid model of compressible particle-laden gas flows is considered. The governing equations are discretized by a high-resolution finite element method based on algebraic flux correction. A multidimensional limiter of TVD type is employed to constrain the local characteristic variables for the continuous gas phase and conservative fluxes for a suspension of solid particles. Special emphasis is laid on the efficient computation of steady state solutions at arbitrary Mach numbers. To avoid stability restrictions and convergence problems, the characteristic boundary conditions are imposed weakly and treated in a fully implicit manner. A two-way coupling via the interphase drag force is implemented using operator splitting. The Douglas-Rachford scheme is found to provide a robust treatment of the interphase exchange terms within the framework of a fractional-step solution strategy. Two-dimensional simulation results are presented for a moving shock wave and for a steady nozzle flow.