Mathematical models and numerical simulations have been prominently successful in providing insight for the exploration of evolution mechanism in scientific engineering. In this paper, a linearized fourth-order compact alternating direction implicit (ADI) method is developed to numerically solve a mathematical model, which includes phytoplankton–zooplankton interactions arising in marine ecosystem. The solvability, convergence and stability of the proposed method are discussed as well. It is proven that the proposed method is fourth-order accurate in space and second-order accurate in time. Finally, we present several numerical simulations to confirm the theoretical results.