This article presents an algorithm for the numerical solution of an initial-boundary value problem for a symmetric t-hyperbolic system of partial differential equations. This problem is based on continual filtration model, which describes the propagation of seismic waves in a poroelastic medium saturated with a fluid characterized by such physical parameters as the propagation velocities of longitudinal P- (fast and slow) and transverse S-waves, the density of the medium materials, and porosity. The system of linearized equations of saturated porous media is formulated in terms of physical variables of the velocity–stress tensor of the porous matrix and the velocity–pressure of the saturating fluid in the absence of energy dissipation. The solution is implemented numerically using an explicit finite difference upwind scheme built on a staggered grid to avoid the appearance of oscillations in the solution functions. The program code implementing parallel computing is developed in the high-performance Julia programming language. The possibility of using the approach is demonstrated by the example of solving the problem of propagation of seismic waves from a source located in the formation. Computational experiments based on real data from oil reservoirs have been implemented, and dynamic visualization of solutions consistent with the first waves arrival times has been obtained.