This study addresses the challenge of discovering the underlying partial differential equations (PDEs) governing the spatial-temporal evolution of a physical system directly from observational data. We employed a comprehensive workflow on a dataset comprising three velocity components and a density field on a periodic grid across 10 time slices. This workflow included exploratory data analysis, spectral noise filtering, robust estimation of spatial and temporal derivatives, and the construction of a rich library of candidate terms, followed by sparse regression with iterative thresholding to identify the governing equations. Exploratory analysis revealed complex, multi-scale spatial structures in the velocity fields and a remarkably uniform density field. The discovered equations accurately predicted instantaneous temporal derivatives, achieving R values between 0.593 and 0.732 for velocity components and 0.362 for density. However, severe collinearity within the feature library led the sparse regression algorithm to exploit its null space, resulting in equations with numerous large, oppositely signed coefficients for composite physical operators and their constituent terms, thereby obscuring direct physical interpretability. Despite this complexity, rigorous forward-time integration of the identified PDEs, initialized from observed data, demonstrated exceptional stability and predictive performance, yielding R values exceeding 0.999 for velocity fields and 0.992 for density over a subsequent time step. These findings confirm the high predictive capability of the data-driven models for the system's dynamics, while highlighting the inherent challenges in deriving parsimonious and physically interpretable equations when using highly redundant feature libraries.