A computational approach with the aid of the Linear Multistep Method (LMM) for the numerical solution of differential equations with initial value problems or boundary conditions has appeared several times in the literature due to its good accuracy and stability properties. The major objective of this article is to extend a multistep approach for the numerical solution of the Partial Differential Equation (PDE) originating from fluid mechanics in a two-dimensional space with initial and boundary conditions, as a result of the importance and utility of the models of partial differential equations in applications, particularly in physical phenomena, such as in convection-diffusion models, and fluid flow problems. Thus, a multistep collocation formula, which is based on orthogonal polynomials is proposed. Ninth-order Multistep Collocation Formulas (NMCFs) were formulated through the principle of interpolation and collocation processes. The theoretical analysis of the NMCFs reveals that they have algebraic order nine, are zero-stable, consistent, and, thus, convergent. The implementation strategies of the NMCFs are comprehensively discussed. Some numerical test problems were presented to evaluate the efficacy and applicability of the proposed formulas. Comparisons with other methods were also presented to demonstrate the new formulas’ productivity. Finally, figures were presented to illustrate the behavior of the numerical examples.