Abstract: In the realm of solving partial differential equations (PDEs), classical numerical methods often require fine mesh grids and reasonable time stepping to meet stability, consistency and convergence conditions. Recently, machine learning has been widely applied to solve PDE problems, but faced with issues of poor interpretability, weak generalizability and strong dependency of rich labeled data. To this end, we introduce a novel PDE-Preserved Coarse Correction Network (P$^2$C$^2$Net) to efficiently solve spatiotemporal PDE systems on coarse mesh grids in small training data regimes. The model consists of two synergistic modules: (1) a trainable PDE block that learns to update the coarse solution, aka., the system state, based on a high-order numerical scheme with boundary condition encoding, and (2) a neural network block that consistently corrects the solution on the fly. In particular, we propose a learnable symmetric Conv filter, with weights shared over the entire model, to accurately estimate the spatial derivatives of PDE based on the neural-corrected system state. The resulting physics-encoded model is capable of handling limited training data (e.g., 3--5 trajectories) and accelerates the prediction of PDE solutions on coarse spatiotemporal grids while maintaining a high accuracy. P$^2$C$^2$Net achieves consistent state-of-the-art performance with over 50\% gain (e.g., in terms of relative prediction error) across four datasets covering complex reaction-diffusion processes and turbulent flows.