Poroacoustic shocks in inhomogeneous gases are numerically simulated using Krylov subspace spectral (KSS) methods; specifically, existing KSS methods are modified to use basis functions tailored to the assumed density profile of the gas, and to handle the nonhomogeneous boundary condition used to insert the shock‐inducing signal. The simulations presented demonstrate that KSS methods are effective at capturing the features of the solution, consistent with theoretical expectations, even when using a time step that greatly exceeds the CFL limit.